{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 高斯玻色采样模拟分子振动激发"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子中振动模式的激发会影响化学反应的结果。\n",
    "振动激发可以为分子提供额外的动能，以克服特定反应的能量壁垒，并有助于控制分子的稳定性。\n",
    "然而，对涉及分子振动和电子态同时变化的过程进行所有振动能级激发概率模拟是具有挑战性的。\n",
    "本例中，我们演示如何使用高斯玻色取样（GBS）来模拟分子振动激发。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 甲酸分子的振动激发"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import deepquantum as dq\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们还需要描述分子振动模式和分子在振动跃迁过程中几何结构变化的分子参数。\n",
    "这些参数由杜辛斯基矩阵和位移向量表示，这些参数是从原子坐标、原子质量、振动频率以及分子的正常模式获得的。\n",
    "这些分子参数可以通过电子结构计算获得。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及$1^2A'$ 态的平衡结构坐标如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ri = np.genfromtxt(\"./data/formic_ri.csv\", delimiter=\",\", skip_header=0)[:, np.newaxis]\n",
    "rf = np.genfromtxt(\"./data/formic_rf.csv\", delimiter=\",\", skip_header=0)[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及 $1^2A'$ 态的质量加权简正坐标如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "li = np.genfromtxt(\"./data/formic_li.csv\", delimiter=\",\", skip_header=0)\n",
    "lf = np.genfromtxt(\"./data/formic_lf.csv\", delimiter=\",\", skip_header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及 $1^2A'$ 态简正模式频率如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = np.genfromtxt(\"./data/formic_omega.csv\", delimiter=\",\", skip_header=0)\n",
    "omegap = np.genfromtxt(\"./data/formic_omegap.csv\", delimiter=\",\", skip_header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算中用到的物理学常量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 299792458.0  # 光速\n",
    "mu = 1.6605390666 * 10**-27  # 原子质量单位\n",
    "h = 6.62607015 * 10**-34  # 普朗克常数\n",
    "\n",
    "m_c = 12  # 碳原子相对原子质量\n",
    "m_h = 1.007825037  # 氢原子相对原子质量\n",
    "m_o = 15.994914640  # 氧原子相对原子质量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Duschinsky 矩阵 $U$ 及位移矢量 $\\delta$ 计算如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.99343181,  0.01440011,  0.01532633,  0.02861045,  0.06378083,\n",
       "         0.07513988, -0.04280796],\n",
       "       [-0.01485231,  0.99314303,  0.07419536,  0.0769291 , -0.03610189,\n",
       "        -0.00248855,  0.01727882],\n",
       "       [-0.01185718, -0.09164895,  0.84227531,  0.17994129, -0.38567948,\n",
       "         0.30738928,  0.08008576],\n",
       "       [ 0.03813492,  0.04087165, -0.34031972, -0.52311845, -0.66785609,\n",
       "         0.38477092,  0.11420873],\n",
       "       [-0.04133251, -0.03419326, -0.40036477,  0.76362651, -0.10356638,\n",
       "         0.48381836,  0.09406589],\n",
       "       [ 0.0907918 , -0.04184572, -0.09066643,  0.31511355, -0.59003401,\n",
       "        -0.719268  ,  0.13037051],\n",
       "       [-0.03245558,  0.00500992, -0.02063661,  0.06935002, -0.20181377,\n",
       "         0.01732396, -0.97588364]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = []\n",
    "for li_ele in li:\n",
    "    for lf_elf in lf:\n",
    "        u.append(np.sum(li_ele * lf_elf))\n",
    "u = np.array(u[-1::-1]).reshape(7, 7).T\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.22536617],\n",
       "       [ 0.14689208],\n",
       "       [ 1.55989779],\n",
       "       [-0.37838396],\n",
       "       [ 0.45525871],\n",
       "       [-0.34391138],\n",
       "       [ 0.06184607]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta = []\n",
    "m = np.diag([m_c, m_c, m_c, m_o, m_o, m_o, m_o, m_o, m_o, m_h, m_h, m_h, m_h, m_h, m_h])\n",
    "for i in range(len(omegap)):\n",
    "    d = lf[i].T @ np.sqrt(m) @ (ri - rf)\n",
    "    l = np.sqrt(h / (4 * np.pi**2 * 100 * omegap[i] * c * mu)) / (10**-10)\n",
    "    delta.append(d / l)\n",
    "delta = np.array(delta[-1::-1])\n",
    "delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Duschinsky矩阵非对角元素描述了不同电子态之间简正模式的混合情况"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAHWCAYAAABKX+K4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0fklEQVR4nO3de3xU1bn/8e9MQiagSbgmBIhEUUDul0gawWpLFNFGqR7LQZQQLiomFcxRMVUB9UioVYS2FASL0P5EaFWs9QLSlKAcoIHQKB6Qm6TkqCRQSgJREsjs3x/K2DSxZGDPrD2Zzzuv/Xoxe/as/czLSx6eZ621XZZlWQIAADDAbToAAAAQvkhEAACAMSQiAADAGBIRAABgDIkIAAAwhkQEAAAYQyICAACMIREBAADGRJoO4Hx4vV599tlniomJkcvlMh0OAKAZsCxLx48fV6dOneR2B/fv6ydPnlRtba1t40VFRSk6Otq28QIhpBORzz77TElJSabDAAA0Q2VlZerSpUvQ7nfy5Em1jLtAqvXaNmbHjh114MABRycjIZ2IxMTEfPWHYQlSZPPuMn3+6jbTIcBGXsu+/9E43SnvKdMhBIWl8PhnGukK6V8bTXL8+HH1uqTvN79jgqS2tvarJGRYRynShir/aUuHNh5SbW0tiUig+Noxke5mn4jExsaaDgE2Cq9ExL4ys5N5wyQRaeFqYTqEoDHW8m9h0+80V2j8OxnSiQgAAM2OW/YsJQmRv5+HSJgAAKA5oiICAICTuFxfHXaMEwJIRAAAcJrQyCFsQWsGAAAYQ0UEAAAnCbPWDBURAABgDBURAACcJMyW75KIAADgJLRmAAAAgoOKCAAATuKSPct3Q6MgQiICAICjuF1fHXaMEwJozQAAAGOoiAAA4CS0ZgAAgDGsmgEAAAgOKiIAADhJmLVmqIgAAABjqIgAAOAkYbZ8l0QEAAAnoTUDAAAQHFREAABwEpbvBt+CBQuUnJys6OhopaamqqioyHRIAACYcWaOiB1HCDCeiKxatUq5ubmaOXOmtm/frv79+2vEiBGqqKgwHRoAAAgw44nI3LlzNXnyZGVlZalXr15atGiRWrVqpaVLl5oODQCA4HPZeIQAo4lIbW2tiouLlZ6e7jvndruVnp6uzZs3N7i+pqZGVVVV9Q4AABC6jCYiR44cUV1dnRISEuqdT0hI0KFDhxpcn5+fr7i4ON+RlJQUrFABAAgOl76ZsHpeh+kv0jTGWzP+yMvLU2Vlpe8oKyszHRIAAPYLk7aMZHj5bvv27RUREaHy8vJ658vLy9WxY8cG13s8Hnk8nmCFBwAAAsxoRSQqKkqDBw9WQUGB75zX61VBQYHS0tIMRgYAgCFhtnzX+IZmubm5yszMVEpKioYMGaJ58+apurpaWVlZpkMDACD4wmyLd+OJyOjRo3X48GHNmDFDhw4d0oABA7RmzZoGE1gBAEDzYzwRkaScnBzl5OSYDgMAAPPY4h0AACA4HFERAQAAX3PLnjJBiJQaSEQAAHASWjMAAADBQUUEAAAnYfkuAAAwhtYMAABAcFARAQDASVg1AwAAjKE1AwAAEBxURAAAcJIwWzVDRQQAABhDRQQAACdxu7467BgnBJCIAADgJExWBQAACA4qIgAAOEmYTVYlEQEAwFFcctnQVrFCJBOhNQMAAIyhIgIAgIO4XPZURORyyTr/UQKOiggAADCGiggAAA5i1+pduRQSFRESEQAAHMRtU2vGcrnktSGeQGsWicjnr25TbGys6TAC6oKRPU2HEBRH3txuOoTgsELh7yn2iIrwmA4hKFq4o0yHEBReKxR+tZ2fSHcL0yGElWaRiAAA0FzYOVk1FJCIAADgIOGWiLBqBgAAGENFBAAAB6EiAgAAECRURAAAcBA79xEJBSQiAAA4CK0ZAACAIKEiAgCAg4RbRYREBAAAB3F9/WPHSKGA1gwAADCGiggAAA5CawYAABgTbst3ac0AAABjqIgAAOAgbpdsac1YVEQAAAD+PRIRAAAc5MxkVTuOc7FgwQIlJycrOjpaqampKioq+rfXz5s3Tz169FDLli2VlJSk+++/XydPnmzy/WjNAADgICZXzaxatUq5ublatGiRUlNTNW/ePI0YMUK7d+9WfHx8g+tXrFihhx9+WEuXLtWVV16pPXv2aPz48XK5XJo7d26T7klFBAAASJLmzp2ryZMnKysrS7169dKiRYvUqlUrLV26tNHrN23apKFDh+r2229XcnKyrrvuOo0ZM+asVZR/RiICAICTuL5Zwns+x5nlu1VVVfWOmpqaRm9bW1ur4uJipaen+8653W6lp6dr8+bNjX7myiuvVHFxsS/x+OSTT/T222/rhhtuaPLXpTUDAICD2NWaOTNGUlJSvfMzZ87UrFmzGlx/5MgR1dXVKSEhod75hIQEffzxx43e4/bbb9eRI0c0bNgwWZal06dP65577tFPfvKTJsdJIgIAQDNWVlam2NhY32uPx2Pb2IWFhZo9e7Z+9atfKTU1Vfv27dPUqVP15JNP6rHHHmvSGCQiAAA4iN0VkdjY2HqJyLdp3769IiIiVF5eXu98eXm5Onbs2OhnHnvsMd15552aNGmSJKlv376qrq7WXXfdpUceeURu99lngDBHBAAAKCoqSoMHD1ZBQYHvnNfrVUFBgdLS0hr9zBdffNEg2YiIiJAkWZbVpPsaTUTee+89ZWRkqFOnTnK5XHr99ddNhgMAgHEu2bSPyDk8bCY3N1dLlizR8uXLtWvXLk2ZMkXV1dXKysqSJI0bN055eXm+6zMyMrRw4UKtXLlSBw4c0Lp16/TYY48pIyPDl5CcjdHWTHV1tfr3768JEybolltuMRkKAACOYHdrxh+jR4/W4cOHNWPGDB06dEgDBgzQmjVrfBNYDx48WK8C8uijj8rlcunRRx/Vp59+qg4dOigjI0NPPfVU0+O0mlo7CTCXy6XVq1dr1KhRTf5MVVWV4uLi9PnfP21S/yuUXTCyp+kQguLIm9tNhxAczvjPLiiiIuybGOdkLdxRpkMICq/lNR1CwFVVVSmxXWdVVlYG9XfLmd9pHfLS5I4+/zqB9+RpHc7fHPTv4a+QmqxaU1NTb/1zVVWVwWgAALCfbx8QG8YJBSE1WTU/P19xcXG+41/XRgMAEOpMP2sm2EIqEcnLy1NlZaXvKCsrMx0SAAA4DyHVmvF4PLZuxAIAgNOYnKxqQkhVRAAAQPNitCJy4sQJ7du3z/f6wIEDKikpUdu2bXXRRRcZjAwAADPcLpfcYTRb1Wgism3bNn3ve9/zvc7NzZUkZWZmatmyZYaiAgDAnHBbNWM0EbnmmmuavAUsAABofkJqsioAAM1duE1WJREBAMBBXDq358Q0Nk4oYNUMAAAwhooIAAAOEm6tGSoiAADAGCoiAAA4SLhVREhEAABwkHDbR4TWDAAAMIaKCAAADkJrBgAAGBNuiQitGQAAYAwVEQAAnMSmikiozFYlEQEAwEFYNQMAABAkVEQAAHAQJqsCAAAECRURAAAc5Ks5InZURGwIJghIRAAAcBBaMwAAAEFCRQQAAAdxyablu+c/RFCQiAAA4CC0ZgAAAIKEikiIqHijyHQIQdH5ketNhxAUHzy23HQIQfM/n28yHUJQJLZKMB1CUFwR/x3TIQTcl6erjd6figgAAECQUBEBAMBBwq0iQiICAICD8NA7AACAIKEiAgCAg9CaAQAA5oRZb4bWDAAAMIaKCAAADhJurRkqIgAAwBgqIgAAOEiYTREhEQEAwElozQAAAAQJFREAABwk3CoiJCIAADhIuCUitGYAAIAxVEQAAHAQVs0AAABjaM0AAAAECRURAACcxKaKSKj0ZqiIAAAAY6iIAADgIOE2R4REBAAABwm3RITWDAAAMIaKCAAADhJu+4gYrYjk5+friiuuUExMjOLj4zVq1Cjt3r3bZEgAABjlksvXnjmvQ6GRiRhNRDZs2KDs7Gxt2bJF69at06lTp3TdddepurraZFgAACBIjLZm1qxZU+/1smXLFB8fr+LiYn33u981FBUAAOaE22RVR80RqayslCS1bdu20fdrampUU1Pje11VVRWUuAAAQGA4ZtWM1+vVtGnTNHToUPXp06fRa/Lz8xUXF+c7kpKSghwlAACBZcv8ELt2Zw0CxyQi2dnZ+uijj7Ry5cpvvSYvL0+VlZW+o6ysLIgRAgAQeGdWzdhxhAJHtGZycnL05ptv6r333lOXLl2+9TqPxyOPxxPEyAAAQCAZTUQsy9KPf/xjrV69WoWFhbr44otNhgMAgHFMVg2i7OxsrVixQn/4wx8UExOjQ4cOSZLi4uLUsmVLk6EBAGCGSzbtaHb+QwSD0TkiCxcuVGVlpa655holJib6jlWrVpkMCwAABInx1gwAAPhGuLVmHLNqBgAAhB9HrJoBAABfcbu+OuwYJxSQiAAA4CC0ZgAAAIKEiggAAA7idrnktqGaYccYwUBFBAAABzH9rJkFCxYoOTlZ0dHRSk1NVVFR0b+9/tixY8rOzlZiYqI8Ho+6d++ut99+u8n3oyICAAAkSatWrVJubq4WLVqk1NRUzZs3TyNGjNDu3bsVHx/f4Pra2lpde+21io+P1yuvvKLOnTvrb3/7m1q3bt3ke5KIAADgIG7Z0644lzHmzp2ryZMnKysrS5K0aNEivfXWW1q6dKkefvjhBtcvXbpUR48e1aZNm9SiRQtJUnJycsDjBAAAAeL6eo7I+R5nWjNVVVX1jpqamkbvW1tbq+LiYqWnp/vOud1upaena/PmzY1+5o033lBaWpqys7OVkJCgPn36aPbs2aqrq2vy9yURAQCgGUtKSlJcXJzvyM/Pb/S6I0eOqK6uTgkJCfXOJyQk+J4F968++eQTvfLKK6qrq9Pbb7+txx57TM8++6z++7//u8nx0ZoBAMBB7N5HpKysTLGxsb7zHo/nvMc+w+v1Kj4+XosXL1ZERIQGDx6sTz/9VD/72c80c+bMJo1BIgIAQDMWGxtbLxH5Nu3bt1dERITKy8vrnS8vL1fHjh0b/UxiYqJatGihiIgI37nLL79chw4dUm1traKios56X79bM+vXr//W955//nl/hwMAAP/Ejvkh57IXSVRUlAYPHqyCggLfOa/Xq4KCAqWlpTX6maFDh2rfvn3yer2+c3v27FFiYmKTkhDpHBKR66+/Xg8++KBOnTrlO3fkyBFlZGQ0OqMWAAA0ncl9RHJzc7VkyRItX75cu3bt0pQpU1RdXe1bRTNu3Djl5eX5rp8yZYqOHj2qqVOnas+ePXrrrbc0e/ZsZWdnN/mefrdm1q9fr3HjxmndunVasWKFDhw4oIkTJ6pHjx4qKSnxdzgAAOAQo0eP1uHDhzVjxgwdOnRIAwYM0Jo1a3wTWA8ePCi3+5saRlJSktauXav7779f/fr1U+fOnTV16lRNnz69yff0OxG58sorVVJSonvuuUeDBg2S1+vVk08+qYceeihkHrADAIBTmdxHRJJycnKUk5PT6HuFhYUNzqWlpWnLli3neLdzjHPPnj3atm2bunTposjISO3evVtffPHFOQcBAAC+YmqOiCl+JyJz5sxRWlqarr32Wn300UcqKirSX//6V/Xr1+9bNzwBAABojN+tmfnz5+v111/XyJEjJUl9+vRRUVGRfvKTn+iaa6751h3bAADA2dm9j4jT+Z2I7NixQ+3bt693rkWLFvrZz36mH/zgB7YFBgAAmj+/E5H27dvr2LFjeuWVV7R//349+OCDatu2rbZv365LL700EDECABA27JrfESpzRPxORD788EOlp6crLi5OpaWlmjx5stq2bavXXntNBw8e1G9+85tAxAkAQFhwfX3YMU4o8DsRuf/++zV+/Hg9/fTTiomJ8Z2/4YYbdPvtt9saXFOdtk7rtHXayL2DxRPZynQIQfH/7n3AdAhBMX/7UtMhBE3x3/7PdAhBcfd30s9+UTNw2nvq7BeFuHD4jk7idyKybds2LV68uMH5zp07f+vT+QAAQNPQmjkLj8ejqqqqBuf37NmjDh062BIUAADhyi2bEpEQac74vY/ITTfdpCeeeML3rBmXy6WDBw9q+vTpuvXWW20PEAAANF9+JyLPPvusTpw4ofj4eH355Ze6+uqrdemllyomJkZPPfVUIGIEACBsmHzonQl+t2bi4uK0bt06bdy4UR9++KFOnDihQYMGKT09PCZqAQAA+/idiJwxbNgwDRs2zM5YAAAIey6bJqs2q4rIz3/+8yYPeN99951zMAAAhDv2EWnEc889V+/14cOH9cUXX6h169aSpGPHjqlVq1aKj48nEQEAAE3WpMmqBw4c8B1PPfWUBgwYoF27duno0aM6evSodu3apUGDBunJJ58MdLwAADRrZ/YRseMIBX6vmnnsscf0i1/8Qj169PCd69Gjh5577jk9+uijtgYHAEC4IRE5i88//1ynTzfcTr2urk7l5eW2BAUAAMKD34nI8OHDdffdd2v79u2+c8XFxZoyZQpLeAEAOE8ul117iZj+Jk3jdyKydOlSdezYUSkpKfJ4PPJ4PBoyZIgSEhL0wgsvBCJGAADCRri1ZvzeR6RDhw56++23tWfPHn388ceSpJ49e6p79+62BwcAAJq3c97QrHv37iQfAADYjH1EzqKurk7Lli1TQUGBKioq5PV6673/5z//2bbgAABA8+Z3IjJ16lQtW7ZMN954o/r06RMyW8gCABAK7Jrf0WzniKxcuVK/+93vdMMNNwQiHgAAwlq4JSJ+r5qJiorSpZdeGohYAABAmPE7Efmv//ovzZ8/X5ZlBSIeAADCmj17iLhCZuqE362ZjRs3av369XrnnXfUu3dvtWjRot77r732mm3BAQAQbtw6hyrBt4wTCvxORFq3bq0f/vCHgYgFAACEGb8TkRdffDEQcQAAAEmyq60SIq2ZUKncAACAZqhJFZFBgwapoKBAbdq00cCBA/9tpvbPD8M7m4ULF2rhwoUqLS2VJPXu3VszZszQyJEjmzwGAADNSbgt321SInLzzTfL4/FIkkaNGmXbzbt06aI5c+bosssuk2VZWr58uW6++Wb99a9/Ve/evW27DwAAoYJEpBEzZ85s9M/nKyMjo97rp556SgsXLtSWLVtIRAAACAPn/NA7u9XV1en3v/+9qqurlZaW1ug1NTU1qqmp8b2uqqoKVngAAASFXXuANNt9ROy2Y8cOpaWl6eTJk7rwwgu1evVq9erVq9Fr8/Pz9fjjjwc5QgAAgsctl9w2PDvXjjGCwfiqmR49eqikpER/+ctfNGXKFGVmZmrnzp2NXpuXl6fKykrfUVZWFuRoAQCAnYxXRP752TWDBw/W1q1bNX/+fD3//PMNrvV4PL5JswAANEfh1po554pIbW2tdu/erdOnT9sZj7xeb715IAAAoPnyOxH54osvNHHiRLVq1Uq9e/fWwYMHJUk//vGPNWfOHL/GysvL03vvvafS0lLt2LFDeXl5Kiws1NixY/0NCwCAZuHM8l07jlDgdyKSl5enDz74QIWFhYqOjvadT09P16pVq/waq6KiQuPGjVOPHj00fPhwbd26VWvXrtW1117rb1gAADQLLht/QoHfc0Ref/11rVq1St/5znfq9Z969+6t/fv3+zXWr3/9a39vDwAAmhG/E5HDhw8rPj6+wfnq6uqQmRgDAIBTMVn1LFJSUvTWW2/5Xp/5oi+88MK3bkQGAACaJtzmiPhdEZk9e7ZGjhypnTt36vTp05o/f7527typTZs2acOGDYGIEQAANFN+V0SGDRumkpISnT59Wn379tW7776r+Ph4bd68WYMHDw5EjAAAhA2Xb2/V8z9CwTltaNatWzctWbLE7lgAAAh7btn09N3mtGrGn4fLxcbGnnMwAAAgvDQpEWndunWTZ9/W1dWdV0AAAIQ1l00rXkKjINK0RGT9+vW+P5eWlurhhx/W+PHjfatkNm/erOXLlys/Pz8wUQIAgGapSYnI1Vdf7fvzE088oblz52rMmDG+czfddJP69u2rxYsXKzMz0/4oAQAIE3btihoqO6v6PaV28+bNSklJaXA+JSVFRUVFtgQFAEC4Crd9RPxORJKSkhpdMfPCCy8oKSnJlqAAAEB48Hv57nPPPadbb71V77zzjlJTUyVJRUVF2rt3r1599VXbAwQAIJywxftZ3HDDDdq7d68yMjJ09OhRHT16VBkZGdqzZ49uuOGGQMQIAEDYcNv4EwrOaUOzLl26aPbs2XbHAgAAwsw5JSLHjh3Tr3/9a+3atUuS1Lt3b02YMEFxcXG2BgcAQLihNXMW27ZtU7du3fTcc8/5WjNz585Vt27dtH379kDECAAAmim/KyL333+/brrpJi1ZskSRkV99/PTp05o0aZKmTZum9957z/YgAQAIF+FWEfE7Edm2bVu9JESSIiMj9dBDDzW6vwgAAGi6M8/OtWOcUOB3ayY2NlYHDx5scL6srEwxMTG2BAUAAMKD34nI6NGjNXHiRK1atUplZWUqKyvTypUrNWnSpHrbvgMAAP+dac3YcYQCv1szzzzzjFwul8aNG6fTp09Lklq0aKEpU6Zozpw5tgcIAEA4sWt79lDZ4t3vRCQqKkrz589Xfn6+9u/fL0nq1q2bWrVqZXtwAACgeTunfUQkqVWrVurbt6+dsZw7y/rqaMZCpcR2vv7+5VHTIQTFm1s+MB1C0Gyd9hvTIQTFvJKFpkMIioiuEaZDCLgIl9nvGG5P321yIjJhwoQmXbd06dJzDgYAAISXJiciy5YtU9euXTVw4EBZzbz6AACAKW6XW27X+T8nxo4xgqHJiciUKVP08ssv68CBA8rKytIdd9yhtm3bBjI2AADCTrhtaNbkdGnBggX6/PPP9dBDD+mPf/yjkpKS9KMf/Uhr166lQgIAAM6JX3Ubj8ejMWPGaN26ddq5c6d69+6te++9V8nJyTpx4kSgYgQAIIy4bPlRc5us+q/cbrdcLpcsy1JdXZ2dMQEAELbCbR8RvyoiNTU1evnll3Xttdeqe/fu2rFjh375y1/q4MGDuvDCCwMVIwAACJIFCxYoOTlZ0dHRSk1NVVFRUZM+t3LlSrlcLo0aNcqv+zW5InLvvfdq5cqVSkpK0oQJE/Tyyy+rffv2ft0MAAD8eyb3EVm1apVyc3O1aNEipaamat68eRoxYoR2796t+Pj4b/1caWmpHnjgAV111VV+37PJiciiRYt00UUX6ZJLLtGGDRu0YcOGRq977bXX/A4CAAB8xe2yp63iPoch5s6dq8mTJysrK0vSV7/733rrLS1dulQPP/xwo5+pq6vT2LFj9fjjj+v999/XsWPH/LpnkxORcePGhcxSIAAA4J/a2loVFxcrLy/Pd87tdis9PV2bN2/+1s898cQTio+P18SJE/X+++/7fV+/NjQDAACB5XK55bJhM7IzY1RVVdU77/F45PF4Glx/5MgR1dXVKSEhod75hIQEffzxx43eY+PGjfr1r3+tkpKSc44zNLZdAwAA5yQpKUlxcXG+Iz8/35Zxjx8/rjvvvFNLliw5rzmj57x8FwAA2M/uyaplZWWKjY31nW+sGiJJ7du3V0REhMrLy+udLy8vV8eOHRtcv3//fpWWliojI8N3zuv1SpIiIyO1e/dudevW7axxkogAAOAgdu8jEhsbWy8R+TZRUVEaPHiwCgoKfEtwvV6vCgoKlJOT0+D6nj17aseOHfXOPfroozp+/Ljmz5+vpKSkJsVJIgIAACRJubm5yszMVEpKioYMGaJ58+apurrat4pm3Lhx6ty5s/Lz8xUdHa0+ffrU+3zr1q0lqcH5f4dEBAAABzH50LvRo0fr8OHDmjFjhg4dOqQBAwZozZo1vgmsBw8elNtt7/RSEhEAABzELZfcNswROdcxcnJyGm3FSFJhYeG//ey5rLBl1QwAADCGiggAAA5isjVjAhURAABgDBURAAAcxO6dVZ2ORAQAAAcxPVk12EIjXQIAAM0SFREAAByEyaqGzJkzRy6XS9OmTTMdCgAABrls+RGtmabbunWrnn/+efXr1890KAAAIIiMJyInTpzQ2LFjtWTJErVp08Z0OAAAGOWSy9eeOa+DikjTZGdn68Ybb1R6evpZr62pqVFVVVW9AwAAhC6jk1VXrlyp7du3a+vWrU26Pj8/X48//niAowIAwByW7wZJWVmZpk6dqpdeeknR0dFN+kxeXp4qKyt9R1lZWYCjBAAguM5saGbHEQqMVUSKi4tVUVGhQYMG+c7V1dXpvffe0y9/+UvV1NQoIiKi3mc8Ho88Hk+wQwUAAAFiLBEZPny4duzYUe9cVlaWevbsqenTpzdIQgAACAffLL89/3FCgbFEJCYmRn369Kl37oILLlC7du0anAcAIFy4XPZsRhYi+5mZXzUDAADCl6O2eC8sLDQdAgAARoVba4aKCAAAMMZRFREAAMJduD30jkQEAAAHYUMzAACAIKEiAgCAg9CaAQAAxri+bs7YMU4oCI0oAQBAs0RFBAAAB6E1AwAAjGFDMwAAgCChIgIAgIO4XS65bWir2DFGMFARAQAAxlARAQDAQcJtjgiJCAAADhJuq2ZozQAAAGOoiAAA4Cj27KwaKrUGEhEAAByE1gwAAECQUBEBAMBB3F8/9s6OcUIBFREAAGAMFREAABwk3OaIkIgAAOAg4bahGa0ZAABgTLOoiNR6a1TjPWk6jIDyuKNNhxAUt3S71XQIQdH59k6mQwiau9bNMh1CUOwp/cx0CEGR3a/GdAgBV+M1+x1pzQAAAGO+asycf8OC1gwAAMBZUBEBAMBB3C6X3Da0VewYIxioiAAAAGOoiAAA4CDhtnyXRAQAAAcJt1UztGYAAIAxVEQAAHAQWjMAAMAYWjMAAABBQkUEAAAHcX/9Y8c4oYBEBAAAB6E1AwAAECRURAAAcJBwWzVDRQQAABhDRQQAACexaY6IQmSOCIkIAAAOQmsGAAAgSKiIAADgIOFWESERAQDASVwue+Z3hMgcEVozAADAGCoiAAA4SLi1ZqiIAAAAY4wmIrNmzfLtqX/m6Nmzp8mQAAAw6l9/L57PEQqMt2Z69+6tP/3pT77XkZHGQwIAwJhwa80Y/60fGRmpjh07mg4DAAAYYHyOyN69e9WpUyddcsklGjt2rA4ePGg6JAAAjHHpm6rI+f2EBqMVkdTUVC1btkw9evTQ559/rscff1xXXXWVPvroI8XExDS4vqamRjU1Nb7XVVVVwQwXAICAc8me+R2hkooYTURGjhzp+3O/fv2Umpqqrl276ne/+50mTpzY4Pr8/Hw9/vjjwQwRAAAEkPHWzD9r3bq1unfvrn379jX6fl5eniorK31HWVlZkCMEACCw7GnLnHtzZsGCBUpOTlZ0dLRSU1NVVFT0rdcuWbJEV111ldq0aaM2bdooPT39317fGEclIidOnND+/fuVmJjY6Psej0exsbH1DgAAYI9Vq1YpNzdXM2fO1Pbt29W/f3+NGDFCFRUVjV5fWFioMWPGaP369dq8ebOSkpJ03XXX6dNPP23yPY0mIg888IA2bNig0tJSbdq0ST/84Q8VERGhMWPGmAwLAABjTFZE5s6dq8mTJysrK0u9evXSokWL1KpVKy1durTR61966SXde++9GjBggHr27KkXXnhBXq9XBQUFTb6n0Tki//d//6cxY8bo73//uzp06KBhw4Zpy5Yt6tChg8mwAAAwxq7NyPwdo7a2VsXFxcrLy/Odc7vdSk9P1+bNm5s0xhdffKFTp06pbdu2Tb6v0URk5cqVJm8PAECz968rTD0ejzweT4Prjhw5orq6OiUkJNQ7n5CQoI8//rhJ95o+fbo6deqk9PT0JsfnqDkiAACEO7tbM0lJSYqLi/Md+fn5AYl7zpw5WrlypVavXq3o6Ogmf874zqoAAOAbdrdmysrK6i3uaKwaIknt27dXRESEysvL650vLy8/6w7ozzzzjObMmaM//elP6tevn19xUhEBAKAZ+9fVpt+WiERFRWnw4MH1JpqemXialpb2reM//fTTevLJJ7VmzRqlpKT4HR8VEQAAHMSuDdrPZYzc3FxlZmYqJSVFQ4YM0bx581RdXa2srCxJ0rhx49S5c2dfe+enP/2pZsyYoRUrVig5OVmHDh2SJF144YW68MILm3RPEhEAABzEZCIyevRoHT58WDNmzNChQ4c0YMAArVmzxjeB9eDBg3K7v2mmLFy4ULW1tfqP//iPeuPMnDlTs2bNatI9SUQAAIBPTk6OcnJyGn2vsLCw3uvS0tLzvh+JCAAADmJqHxFTmKwKAACMoSICAICDmJwjYgKJCAAADhJuiQitGQAAYAwVEQAAnMSmyaoKkcmqJCIAADiK6+vDjnGcj9YMAAAwhooIAAAOwj4iAAAAQUJFBAAABwm35bskIgAAOEi4JSK0ZgAAgDFURAAAcJBwm6xKIgIAgIN8tYuIHa2Z0EBrBgAAGENFBAAABwm3yarNIhGJjmiplhGtTIcRUJHuFqZDCIoId7P4V/Ksruw4zHQIQZMc09V0CEFxyQ09TIcQFBfc1Mt0CIF3yms6grASHv/XBwAgRDBZFQAAGBNurRkmqwIAAGOoiAAA4CC0ZgAAgDG0ZgAAAIKEiggAAI7ikj37ooZGRYREBAAABwmvNITWDAAAMIiKCAAADhJuq2aoiAAAAGOoiAAA4CjhNUuERAQAAAcJrzSE1gwAADCIiggAAI4SXjUREhEAAByEVTMAAABBQiICAACMIREBAADGMEcEAAAHcX39Y8c4oYBEBAAABwm3RITWDAAAMIZEBAAAGGM8Efn00091xx13qF27dmrZsqX69u2rbdu2mQ4LAAAjzuwjYscRCozOEfnHP/6hoUOH6nvf+57eeecddejQQXv37lWbNm1MhgUAAILEaCLy05/+VElJSXrxxRd95y6++GKDEQEAgGAy2pp54403lJKSottuu03x8fEaOHCglixZYjIkAAAQREYTkU8++UQLFy7UZZddprVr12rKlCm67777tHz58kavr6mpUVVVVb0DAIDmxWXLDw+9awKv16uUlBTNnj1bkjRw4EB99NFHWrRokTIzMxtcn5+fr8cffzzYYQIAEETh9fRdoxWRxMRE9erVq965yy+/XAcPHmz0+ry8PFVWVvqOsrKyYIQJAAACxGhFZOjQodq9e3e9c3v27FHXrl0bvd7j8cjj8QQjNAAAjAiveojhROT+++/XlVdeqdmzZ+tHP/qRioqKtHjxYi1evNhkWAAAGGPXHiChso+I0dbMFVdcodWrV+vll19Wnz599OSTT2revHkaO3asybAAAECQGH/o3Q9+8AP94Ac/MB0GAAAOEV7NGeOJCAAA+EZ4pSEOeNYMAAAIX1REAABwnFCpZ5w/KiIAAMAYKiIAADgIy3cBAACChEQEAAAYQ2sGAAAH+ebpuec/TiggEQEAwFHCaycRWjMAAMAYKiIAADhIeNVDqIgAAACDqIgAAOAg4baPCIkIAACOEl7NGVozAADAGCoiAAA4SHjVQ0hEAABwmPBKRWjNAAAAY6iIAADgIOG2aoaKCAAAMIZEBAAAGENrBgAAB+HpuyHEsixJ0vHjxw1HEniR7hamQwiKOqvOdAhBUVN30nQIQXOi+oTpEIKiSlWmQwiOU17TEQTe6a++45nfMcFWVWXP7zS7xgm0kE5EziQgPS/ubTgSAEBzc/z4ccXFxQXtflFRUerYsaMuS+5u25gdO3ZUVFSUbeMFgssylfLZwOv16rPPPlNMTEzQZgdXVVUpKSlJZWVlio2NDco9TeB7Nj/h8l35ns2Lie9pWZaOHz+uTp06ye0O7lTKkydPqra21rbxoqKiFB0dbdt4gRDSFRG3260uXboYuXdsbGyz/o//DL5n8xMu35Xv2bwE+3sGsxLyz6Kjox2fONiNVTMAAMAYEhEAAGAMiYifPB6PZs6cKY/HYzqUgOJ7Nj/h8l35ns1LuHzPcBbSk1UBAEBooyICAACMIREBAADGkIgAAABjSET8tGDBAiUnJys6OlqpqakqKioyHZKt3nvvPWVkZKhTp05yuVx6/fXXTYcUEPn5+briiisUExOj+Ph4jRo1Srt37zYdlu0WLlyofv36+fZgSEtL0zvvvGM6rICbM2eOXC6Xpk2bZjoUW82aNcv3iPgzR8+ePU2HFRCffvqp7rjjDrVr104tW7ZU3759tW3bNtNhIQBIRPywatUq5ebmaubMmdq+fbv69++vESNGqKKiwnRotqmurlb//v21YMEC06EE1IYNG5Sdna0tW7Zo3bp1OnXqlK677jpVV1ebDs1WXbp00Zw5c1RcXKxt27bp+9//vm6++Wb97//+r+nQAmbr1q16/vnn1a9fP9OhBETv3r31+eef+46NGzeaDsl2//jHPzR06FC1aNFC77zzjnbu3Klnn31Wbdq0MR0aAsFCkw0ZMsTKzs72va6rq7M6depk5efnG4wqcCRZq1evNh1GUFRUVFiSrA0bNpgOJeDatGljvfDCC6bDCIjjx49bl112mbVu3Trr6quvtqZOnWo6JFvNnDnT6t+/v+kwAm769OnWsGHDTIeBIKEi0kS1tbUqLi5Wenq675zb7VZ6ero2b95sMDLYobKyUpLUtm1bw5EETl1dnVauXKnq6mqlpaWZDicgsrOzdeONN9b777S52bt3rzp16qRLLrlEY8eO1cGDB02HZLs33nhDKSkpuu222xQfH6+BAwdqyZIlpsNCgJCINNGRI0dUV1enhISEeucTEhJ06NAhQ1HBDl6vV9OmTdPQoUPVp08f0+HYbseOHbrwwgvl8Xh0zz33aPXq1erVq5fpsGy3cuVKbd++Xfn5+aZDCZjU1FQtW7ZMa9as0cKFC3XgwAFdddVVvieRNxeffPKJFi5cqMsuu0xr167VlClTdN9992n58uWmQ0MAhPRD7wA7ZGdn66OPPmqWvXZJ6tGjh0pKSlRZWalXXnlFmZmZ2rBhQ7NKRsrKyjR16lStW7euWT8wbOTIkb4/9+vXT6mpqeratat+97vfaeLEiQYjs5fX61VKSopmz54tSRo4cKA++ugjLVq0SJmZmYajg92oiDRR+/btFRERofLy8nrny8vL1bFjR0NR4Xzl5OTozTff1Pr16409yTnQoqKidOmll2rw4MHKz89X//79NX/+fNNh2aq4uFgVFRUaNGiQIiMjFRkZqQ0bNujnP/+5IiMjVVdXZzrEgGjdurW6d++uffv2mQ7FVomJiQ0S5csvv7xZtqFAItJkUVFRGjx4sAoKCnznvF6vCgoKmm2/vTmzLEs5OTlavXq1/vznP+viiy82HVLQeL1e1dTUmA7DVsOHD9eOHTtUUlLiO1JSUjR27FiVlJQoIiLCdIgBceLECe3fv1+JiYmmQ7HV0KFDGyyn37Nnj7p27WooIgQSrRk/5ObmKjMzUykpKRoyZIjmzZun6upqZWVlmQ7NNidOnKj3t6sDBw6opKREbdu21UUXXWQwMntlZ2drxYoV+sMf/qCYmBjfPJ+4uDi1bNnScHT2ycvL08iRI3XRRRfp+PHjWrFihQoLC7V27VrTodkqJiamwfyeCy64QO3atWtW834eeOABZWRkqGvXrvrss880c+ZMRUREaMyYMaZDs9X999+vK6+8UrNnz9aPfvQjFRUVafHixVq8eLHp0BAIppfthJpf/OIX1kUXXWRFRUVZQ4YMsbZs2WI6JFutX7/ektTgyMzMNB2arRr7jpKsF1980XRotpowYYLVtWtXKyoqyurQoYM1fPhw69133zUdVlA0x+W7o0ePthITE62oqCirc+fO1ujRo619+/aZDisg/vjHP1p9+vSxPB6P1bNnT2vx4sWmQ0KA8PRdAABgDHNEAACAMSQiAADAGBIRAABgDIkIAAAwhkQEAAAYQyICAACMIREBAADGkIgAAABjSESAEFNYWCiXy6Vjx46d8xizZs3SgAEDzjuW5ORkzZs377zHARC+SEQAG40fP14ul0v33HNPg/eys7Plcrk0fvz44Af2Lx544IF6D3AEAFNIRACbJSUlaeXKlfryyy99506ePKkVK1Y45sGBF154odq1a2c6DAAgEQHsNmjQICUlJem1117znXvttdd00UUXaeDAgfWuramp0X333af4+HhFR0dr2LBh2rp1a71r3n77bXXv3l0tW7bU9773PZWWlja458aNG3XVVVepZcuWSkpK0n333afq6upvjfFfWzPjx4/XqFGj9MwzzygxMVHt2rVTdna2Tp065bumoqJCGRkZatmypS6++GK99NJLDcY9duyYJk2apA4dOig2Nlbf//739cEHH0iSDh8+rI4dO2r27Nm+6zdt2qSoqCiqM0AYIxEBAmDChAl68cUXfa+XLl2qrKysBtc99NBDevXVV7V8+XJt375dl156qUaMGKGjR49KksrKynTLLbcoIyNDJSUlmjRpkh5++OF6Y+zfv1/XX3+9br31Vn344YdatWqVNm7cqJycHL9iXr9+vfbv36/169dr+fLlWrZsmZYtW+Z7f/z48SorK9P69ev1yiuv6Fe/+pUqKirqjXHbbbepoqJC77zzjoqLizVo0CANHz5cR48eVYcOHbR06VLNmjVL27Zt0/Hjx3XnnXcqJydHw4cP9ytWAM2I6cf/As1JZmamdfPNN1sVFRWWx+OxSktLrdLSUis6Oto6fPiwdfPNN1uZmZmWZVnWiRMnrBYtWlgvvfSS7/O1tbVWp06drKefftqyLMvKy8uzevXqVe8e06dPtyRZ//jHPyzLsqyJEydad911V71r3n//fcvtdltffvllo3HOnDnT6t+/f724u3btap0+fdp37rbbbrNGjx5tWZZl7d6925JkFRUV+d7ftWuXJcl67rnnfPeMjY21Tp48We9e3bp1s55//nnf63vvvdfq3r27dfvtt1t9+/ZtcD2A8BJpOA8CmqUOHTroxhtv1LJly2RZlm688Ua1b9++3jX79+/XqVOnNHToUN+5Fi1aaMiQIdq1a5ckadeuXUpNTa33ubS0tHqvP/jgA3344Yf1WiWWZcnr9erAgQO6/PLLmxRz7969FRER4XudmJioHTt2+OKIjIzU4MGDfe/37NlTrVu3rhfHiRMnGsw9+fLLL7V//37f62eeeUZ9+vTR73//exUXF8vj8TQpPgDNE4kIECATJkzwtUcWLFgQsPucOHFCd999t+67774G7/kzObZFixb1XrtcLnm9Xr/iSExMVGFhYYP3/jlh2b9/vz777DN5vV6Vlpaqb9++Tb4HgOaHRAQIkOuvv161tbVyuVwaMWJEg/e7deumqKgo/c///I+6du0qSTp16pS2bt2qadOmSZIuv/xyvfHGG/U+t2XLlnqvBw0apJ07d+rSSy8NzBfRV9WP06dPq7i4WFdccYUkaffu3fX2Mhk0aJAOHTqkyMhIJScnNzpObW2t7rjjDo0ePVo9evTQpEmTtGPHDsXHxwcsdgDOxmRVIEAiIiK0a9cu7dy5s17L44wLLrhAU6ZM0YMPPqg1a9Zo586dmjx5sr744gtNnDhRknTPPfdo7969evDBB7V7926tWLGi3gRSSZo+fbo2bdqknJwclZSUaO/evfrDH/7g92TVf6dHjx66/vrrdffdd+svf/mLiouLNWnSJLVs2dJ3TXp6utLS0jRq1Ci9++67Ki0t1aZNm/TII49o27ZtkqRHHnlElZWV+vnPf67p06ere/fumjBhgm1xAgg9JCJAAMXGxio2NvZb358zZ45uvfVW3XnnnRo0aJD27duntWvXqk2bNpK+aq28+uqrev3119W/f38tWrSo3vJXSerXr582bNigPXv26KqrrtLAgQM1Y8YMderUydbv8uKLL6pTp066+uqrdcstt+iuu+6qV8lwuVx6++239d3vfldZWVnq3r27/vM//1N/+9vflJCQoMLCQs2bN0+//e1vFRsbK7fbrd/+9rd6//33tXDhQltjBRA6XJZlWaaDAAAA4YmKCAAAMIZEBAAAGEMiAgAAjCERAQAAxpCIAAAAY0hEAACAMSQiAADAGBIRAABgDIkIAAAwhkQEAAAYQyICAACMIREBAADG/H+u0CPfPWuRRwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(abs(u), cmap=\"Greens\")\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"Mode index\")\n",
    "plt.ylabel(\"Mode index\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "非对角元素的存在显示了不同电子态之间简正模式的混合。\n",
    "由 Duschinsky 矩阵和位移向量计算 GBS 参数如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre_transition_squeezing = np.sqrt(omega[-1::-1])\n",
    "post_transition_squeezing = np.sqrt(omegap[-1::-1])\n",
    "\n",
    "j_mat = (\n",
    "    np.diag(post_transition_squeezing)\n",
    "    @ u\n",
    "    @ np.linalg.inv(np.diag(pre_transition_squeezing))\n",
    ")\n",
    "\n",
    "cl, lambda_1, cr = np.linalg.svd(j_mat)\n",
    "\n",
    "delta_2 = np.linalg.inv(j_mat) @ delta / np.sqrt(2)\n",
    "delta_2 = delta_2.flatten()\n",
    "lambda_2 = np.log(lambda_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们现在可以计算每个振动模式中的平均振动量子数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "modes = 7  # 简正模式数量\n",
    "cutoff = 3\n",
    "shots = 500000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"6.363636363636363cm\" version=\"1.1\" width=\"12.0cm\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><polyline fill=\"none\" points=\"40,30 130,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"25\" /><text font-size=\"9\" x=\"83\" y=\"20\">D</text><text font-size=\"7\" x=\"95\" y=\"18\">r =0.056</text><text font-size=\"7\" x=\"95\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"40,60 130,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"55\" /><text font-size=\"9\" x=\"83\" y=\"50\">D</text><text font-size=\"7\" x=\"95\" y=\"48\">r =-0.053</text><text font-size=\"7\" x=\"95\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"40,90 130,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"85\" /><text font-size=\"9\" x=\"83\" y=\"80\">D</text><text font-size=\"7\" x=\"95\" y=\"78\">r =0.982</text><text font-size=\"7\" x=\"95\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"40,120 130,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"115\" /><text font-size=\"9\" x=\"83\" y=\"110\">D</text><text font-size=\"7\" x=\"95\" y=\"108\">r =0.525</text><text font-size=\"7\" x=\"95\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"40,150 130,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"145\" /><text font-size=\"9\" x=\"83\" y=\"140\">D</text><text font-size=\"7\" x=\"95\" y=\"138\">r =-0.112</text><text font-size=\"7\" x=\"95\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"40,180 130,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"175\" /><text font-size=\"9\" x=\"83\" y=\"170\">D</text><text font-size=\"7\" x=\"95\" y=\"168\">r =0.525</text><text font-size=\"7\" x=\"95\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"40,210 130,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"205\" /><text font-size=\"9\" x=\"83\" y=\"200\">D</text><text font-size=\"7\" x=\"95\" y=\"198\">r =-0.016</text><text font-size=\"7\" x=\"95\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"130,30 150,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,30 220,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,60 150,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,60 220,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,90 150,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,90 220,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,120 150,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,120 220,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,150 150,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,150 220,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,180 150,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,180 220,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,210 150,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,210 220,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"150\" y=\"20\" /><text font-size=\"10\" x=\"171\" y=\"120.0\">U</text><polyline fill=\"none\" points=\"220,30 310,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"25\" /><text font-size=\"9\" x=\"263\" y=\"20\">S</text><text font-size=\"7\" x=\"275\" y=\"18\">r =-0.097</text><text font-size=\"7\" x=\"275\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"220,60 310,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"55\" /><text font-size=\"9\" x=\"263\" y=\"50\">S</text><text font-size=\"7\" x=\"275\" y=\"48\">r =-0.07</text><text font-size=\"7\" x=\"275\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"220,90 310,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"85\" /><text font-size=\"9\" x=\"263\" y=\"80\">S</text><text font-size=\"7\" x=\"275\" y=\"78\">r =-0.021</text><text font-size=\"7\" x=\"275\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"220,120 310,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"115\" /><text font-size=\"9\" x=\"263\" y=\"110\">S</text><text font-size=\"7\" x=\"275\" y=\"108\">r =0.06</text><text font-size=\"7\" x=\"275\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"220,150 310,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"145\" /><text font-size=\"9\" x=\"263\" y=\"140\">S</text><text font-size=\"7\" x=\"275\" y=\"138\">r =0.075</text><text font-size=\"7\" x=\"275\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"220,180 310,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"175\" /><text font-size=\"9\" x=\"263\" y=\"170\">S</text><text font-size=\"7\" x=\"275\" y=\"168\">r =0.112</text><text font-size=\"7\" x=\"275\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"220,210 310,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"205\" /><text font-size=\"9\" x=\"263\" y=\"200\">S</text><text font-size=\"7\" x=\"275\" y=\"198\">r =0.187</text><text font-size=\"7\" x=\"275\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"310,30 330,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,30 400,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,60 330,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,60 400,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,90 330,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,90 400,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,120 330,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,120 400,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,150 330,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,150 400,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,180 330,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,180 400,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,210 330,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,210 400,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"330\" y=\"20\" /><text font-size=\"10\" x=\"351\" y=\"120.0\">U</text><text font-size=\"12\" x=\"25\" y=\"30\">0</text><text font-size=\"12\" x=\"25\" y=\"60\">1</text><text font-size=\"12\" x=\"25\" y=\"90\">2</text><text font-size=\"12\" x=\"25\" y=\"120\">3</text><text font-size=\"12\" x=\"25\" y=\"150\">4</text><text font-size=\"12\" x=\"25\" y=\"180\">5</text><text font-size=\"12\" x=\"25\" y=\"210\">6</text></svg>"
      ],
      "text/plain": [
       "<svgwrite.drawing.Drawing at 0x155587650>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cir = dq.photonic.QumodeCircuit(\n",
    "    nmode=modes,\n",
    "    init_state=\"vac\",\n",
    "    # init_state=init_state,\n",
    "    cutoff=cutoff,\n",
    "    backend=\"gaussian\",\n",
    ")\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.d(wires=[i], r=delta_2[i])\n",
    "\n",
    "cir.any(cr, wires=list(range(modes)))\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.s(wires=[i], r=-lambda_2[i])\n",
    "\n",
    "cir.any(cl, wires=list(range(modes)))\n",
    "\n",
    "state = cir()\n",
    "\n",
    "# 线路可视化\n",
    "cir.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n",
      "chain 1: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:03<00:00, 26294.68it/s]\u001b[0m\n",
      "chain 2: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 46735.47it/s]\u001b[0m\n",
      "chain 3: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:01<00:00, 50268.94it/s]\u001b[0m\n",
      "chain 4: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 49138.20it/s]\u001b[0m\n",
      "chain 5: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 49228.56it/s]\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sample = cir.measure(shots=shots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGGCAYAAADmRxfNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABH1UlEQVR4nO3dd3hUZf7//9dMKpBCT0JvLkVKIEKIsi4gJGosrMBXRCmCuCA9GjCuNFe/IK4KiojrKuBnFwU+WChLMxRFIkgwdFABKabAAklIMIEk9+8Pf8yXEZCMk+SQyfNxXXNdnnPfc87LWxzmPeec+7YZY4wAAAAAwA12qwMAAAAAKP8oLAAAAAC4jcICAAAAgNsoLAAAAAC4jcICAAAAgNsoLAAAAAC4jcICAAAAgNsoLAAAAAC4zdvqABVJUVGRUlNTFRgYKJvNZnUcAAAA4DcZY3T+/HnVqVNHdvtvX5OgsChDqampql+/vtUxAAAAAJecOHFC9erV+80+FBZlKDAwUNIv/2GCgoIsTgMAAAD8tuzsbNWvX9/xPfa3UFiUocu3PwUFBVFYAAAAoNwozm38PLwNAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG3eVgcAgJJgm2azOkK5YKYYqyMAADwUVywAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuO2mKSxmzJghm82mcePGOfbl5eVp5MiRqlGjhgICAtS7d29lZGQ4ve/48eOKjY1V5cqVVbt2bcXHx6ugoMCpz6ZNm9ShQwf5+fmpWbNmWrBgwVXnf+utt9SoUSP5+/srMjJS27dvd2ovThYAAACgoropCotvvvlG77zzjtq2beu0f/z48VqxYoWWLl2qzZs3KzU1VQ899JCjvbCwULGxsbp48aK2bt2qhQsXasGCBZo8ebKjz9GjRxUbG6tu3bopJSVF48aN0xNPPKG1a9c6+ixevFhxcXGaMmWKdu7cqXbt2ikmJkanTp0qdhYAAACgIrMZY4yVAXJyctShQwfNnTtXL774osLDwzVr1ixlZWWpVq1aWrRokfr06SNJOnjwoFq2bKmkpCR17txZq1ev1n333afU1FSFhIRIkubNm6eJEyfq9OnT8vX11cSJE7Vq1Srt3bvXcc5+/fopMzNTa9askSRFRkaqY8eOmjNnjiSpqKhI9evX1+jRo/Xss88WK0txZGdnKzg4WFlZWQoKCiqxMQQg2abZrI5QLpgpln7kAwDKGVe+v1p+xWLkyJGKjY1Vjx49nPYnJyfr0qVLTvtbtGihBg0aKCkpSZKUlJSkNm3aOIoKSYqJiVF2drb27dvn6PPrY8fExDiOcfHiRSUnJzv1sdvt6tGjh6NPcbIAAAAAFZm3lSf/6KOPtHPnTn3zzTdXtaWnp8vX11dVq1Z12h8SEqL09HRHnyuLisvtl9t+q092drZ+/vlnnTt3ToWFhdfsc/DgwWJnuZb8/Hzl5+c7trOzs6/bFwAAACjPLLticeLECY0dO1b//ve/5e/vb1WMUjV9+nQFBwc7XvXr17c6EgAAAFAqLCsskpOTderUKXXo0EHe3t7y9vbW5s2b9cYbb8jb21shISG6ePGiMjMznd6XkZGh0NBQSVJoaOhVMzNd3r5Rn6CgIFWqVEk1a9aUl5fXNftceYwbZbmWhIQEZWVlOV4nTpwo3uAAAAAA5YxlhcVdd92lPXv2KCUlxfG67bbb9Oijjzr+2cfHR4mJiY73HDp0SMePH1dUVJQkKSoqSnv27HGavWn9+vUKCgpSq1atHH2uPMblPpeP4evrq4iICKc+RUVFSkxMdPSJiIi4YZZr8fPzU1BQkNMLAAAA8ESWPWMRGBio1q1bO+2rUqWKatSo4dg/dOhQxcXFqXr16goKCtLo0aMVFRXlmIUpOjparVq10oABAzRz5kylp6fr+eef18iRI+Xn5ydJGj58uObMmaMJEyZoyJAh2rBhg5YsWaJVq1Y5zhsXF6dBgwbptttuU6dOnTRr1izl5ubq8ccflyQFBwffMAsAAABQkVn68PaNvP7667Lb7erdu7fy8/MVExOjuXPnOtq9vLy0cuVKjRgxQlFRUapSpYoGDRqkF154wdGncePGWrVqlcaPH6/Zs2erXr16+uc//6mYmBhHn4cfflinT5/W5MmTlZ6ervDwcK1Zs8bpge4bZQEAAAAqMsvXsahIWMcCKD2sY1E8rGMBAHBFuVrHAgAAAED5R2EBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADcRmEBAAAAwG0UFgAAAADc5lJhcenSJTVt2lQHDhworTwAAAAAyiGXCgsfHx/l5eWVVhYAAAAA5ZTLt0KNHDlSL7/8sgoKCkojDwAAAIByyNvVN3zzzTdKTEzUunXr1KZNG1WpUsWp/eOPPy6xcAAAAADKB5cLi6pVq6p3796lkQUAAABAOeVyYTF//vzSyAEAAACgHPtd080WFBTo888/1zvvvKPz589LklJTU5WTk1Oi4QAAAACUDy5fsTh27JjuvvtuHT9+XPn5+erZs6cCAwP18ssvKz8/X/PmzSuNnAAAAABuYi5fsRg7dqxuu+02nTt3TpUqVXLs//Of/6zExMQSDQcAAACgfHD5isWXX36prVu3ytfX12l/o0aN9NNPP5VYMAAAAADlh8tXLIqKilRYWHjV/pMnTyowMLBEQgEAAAAoX1wuLKKjozVr1izHts1mU05OjqZMmaJ77723JLMBAAAAKCdcvhXq1VdfVUxMjFq1aqW8vDz1799f33//vWrWrKkPP/ywNDICAAAAuMm5XFjUq1dPu3bt0kcffaTdu3crJydHQ4cO1aOPPur0MDcAAACAisPlwkKSvL299dhjj5V0FgAAAADl1O8qLA4dOqQ333xTBw4ckCS1bNlSo0aNUosWLUo0HAAAAIDyweWHt5ctW6bWrVsrOTlZ7dq1U7t27bRz5061adNGy5YtK42MAAAAAG5yLl+xmDBhghISEvTCCy847Z8yZYomTJig3r17l1g4AAAAAOWDy1cs0tLSNHDgwKv2P/bYY0pLS3PpWG+//bbatm2roKAgBQUFKSoqSqtXr3a05+XlaeTIkapRo4YCAgLUu3dvZWRkOB3j+PHjio2NVeXKlVW7dm3Fx8eroKDAqc+mTZvUoUMH+fn5qVmzZlqwYMFVWd566y01atRI/v7+ioyM1Pbt253ai5MFAAAAqKhcLiy6du2qL7/88qr9W7Zs0R//+EeXjlWvXj3NmDFDycnJ2rFjh7p3764HH3xQ+/btkySNHz9eK1as0NKlS7V582alpqbqoYcecry/sLBQsbGxunjxorZu3aqFCxdqwYIFmjx5sqPP0aNHFRsbq27duiklJUXjxo3TE088obVr1zr6LF68WHFxcZoyZYp27typdu3aKSYmRqdOnXL0uVEWAAAAoCKzGWPMjTotX77c8c+pqamaPHmy/s//+T/q3LmzJOnrr7/W0qVLNW3aNA0fPtytQNWrV9crr7yiPn36qFatWlq0aJH69OkjSTp48KBatmyppKQkde7cWatXr9Z9992n1NRUhYSESJLmzZuniRMn6vTp0/L19dXEiRO1atUq7d2713GOfv36KTMzU2vWrJEkRUZGqmPHjpozZ46kX1YXr1+/vkaPHq1nn31WWVlZN8xSHNnZ2QoODlZWVpaCgoLcGicAzmzTbFZHKBfMlBt+5AMA4ODK99diPWPRq1evq/bNnTtXc+fOddo3cuTI311YFBYWaunSpcrNzVVUVJSSk5N16dIl9ejRw9GnRYsWatCggePLfFJSktq0aeMoKiQpJiZGI0aM0L59+9S+fXslJSU5HeNyn3HjxkmSLl68qOTkZCUkJDja7Xa7evTooaSkJEkqVpZryc/PV35+vmM7Ozv7d40NAAAAcLMr1q1QRUVFxXoVFha6HGDPnj0KCAiQn5+fhg8frk8++UStWrVSenq6fH19VbVqVaf+ISEhSk9PlySlp6c7FRWX2y+3/Vaf7Oxs/fzzz/rvf/+rwsLCa/a58hg3ynIt06dPV3BwsONVv3794g0KAAAAUM64/IxFSWvevLlSUlK0bds2jRgxQoMGDdL+/futjlUiEhISlJWV5XidOHHC6kgAAABAqfhdC+R988032rhxo06dOqWioiKnttdee82lY/n6+qpZs2aSpIiICH3zzTeaPXu2Hn74YV28eFGZmZlOVwoyMjIUGhoqSQoNDb1q9qbLMzVd2efXszdlZGQoKChIlSpVkpeXl7y8vK7Z58pj3CjLtfj5+cnPz8+F0QAAAADKJ5evWPzf//t/FRkZqfnz52vHjh369ttvHa+UlBS3AxUVFSk/P18RERHy8fFRYmKio+3QoUM6fvy4oqKiJElRUVHas2eP0+xN69evV1BQkFq1auXoc+UxLve5fAxfX19FREQ49SkqKlJiYqKjT3GyAAAAABWZy1csZs+erffff1+DBw92++QJCQm655571KBBA50/f16LFi3Spk2btHbtWgUHB2vo0KGKi4tT9erVFRQUpNGjRysqKsrxsHR0dLRatWqlAQMGaObMmUpPT9fzzz+vkSNHOq4UDB8+XHPmzNGECRM0ZMgQbdiwQUuWLNGqVascOeLi4jRo0CDddttt6tSpk2bNmqXc3Fw9/vjjklSsLAAAAEBF5nJhYbfbdccdd5TIyU+dOqWBAwcqLS1NwcHBatu2rdauXauePXtKkl5//XXZ7Xb17t1b+fn5iomJcZqJysvLSytXrtSIESMUFRWlKlWqaNCgQU6rgjdu3FirVq3S+PHjNXv2bNWrV0///Oc/FRMT4+jz8MMP6/Tp05o8ebLS09MVHh6uNWvWOD3QfaMsAAAAQEVWrHUsrjRz5kylpqZq1qxZpRTJc7GOBVB6WMeieFjHAgDgihJfx+JKzzzzjGJjY9W0aVO1atVKPj4+Tu0ff/yxq4cEAAAAUM65XFiMGTNGGzduVLdu3VSjRg3ZbPxKCAAAAFR0LhcWCxcu1LJlyxQbG1saeQAAAACUQy5PN1u9enU1bdq0NLIAAAAAKKdcLiymTp2qKVOm6MKFC6WRBwAAAEA55PKtUG+88YYOHz6skJAQNWrU6KqHt3fu3Fli4QAAAACUDy4XFr169SqFGAAAAADKM5cLiylTppRGDgAAAADlmMvPWAAAAADAr7l8xcJut//m2hWFhYVuBQIAAABQ/rhcWHzyySdO25cuXdK3336rhQsXatq0aSUWDAAAAED54XJh8eCDD161r0+fPrr11lu1ePFiDR06tESCAQAAACg/SuwZi86dOysxMbGkDgcAAACgHCmRwuLnn3/WG2+8obp165bE4QAAAACUMy7fClWtWjWnh7eNMTp//rwqV66sf/3rXyUaDgAAAED54HJhMWvWLKdtu92uWrVqKTIyUtWqVSupXAAAAADKEZcLi0GDBpVGDgAAAADlmMuFhSRlZmZq+/btOnXqlIqKipzaBg4cWCLBAAAAAJQfLhcWK1as0KOPPqqcnBwFBQU5PW9hs9koLAAAAIAKyOVZoZ5++mkNGTJEOTk5yszM1Llz5xyvs2fPlkZGAAAAADc5lwuLn376SWPGjFHlypVLIw8AAACAcsjlwiImJkY7duwojSwAAAAAyimXn7GIjY1VfHy89u/frzZt2sjHx8ep/YEHHiixcAAAAADKB5sxxrjyBrv9+hc5bDabCgsL3Q7lqbKzsxUcHKysrCwFBQVZHQfwKLZptht3gswUlz7yAQAVnCvfX12+YvHr6WUBAAAAwOVnLAAAAADg1ygsAAAAALiNwgIAAACA2ygsAAAAALiNwgIAAACA21yeFUr6ZWaoH374QadOnbpqlqg777yzRIIBAAAAKD9cLiy+/vpr9e/fX8eOHdOvl8BgHQsAAACgYnK5sBg+fLhuu+02rVq1SmFhYbLZWJQKAAAAqOhcLiy+//57/e///q+aNWtWGnkAAAAAlEMuP7wdGRmpH374oTSyAAAAACinXL5iMXr0aD399NNKT09XmzZt5OPj49Tetm3bEgsHAAAAoHxwubDo3bu3JGnIkCGOfTabTcYYHt4GAAAAKiiXC4ujR4+WRg4AAAAA5ZjLhUXDhg1LIwcAAACAcux3LZB3+PBhzZo1SwcOHJAktWrVSmPHjlXTpk1LNBwAAACA8sHlWaHWrl2rVq1aafv27Wrbtq3atm2rbdu26dZbb9X69etLIyMAAACAm5zLVyyeffZZjR8/XjNmzLhq/8SJE9WzZ88SCwcAAACgfHD5isWBAwc0dOjQq/YPGTJE+/fvL5FQAAAAAMoXlwuLWrVqKSUl5ar9KSkpql27dklkAgAAAFDOuHwr1LBhw/Tkk0/qyJEjuv322yVJX331lV5++WXFxcWVeEAAAAAANz+XC4tJkyYpMDBQr776qhISEiRJderU0dSpUzVmzJgSDwgAAADg5ufyrVA2m03jx4/XyZMnlZWVpaysLJ08eVJjx46VzWZz6VjTp09Xx44dFRgYqNq1a6tXr146dOiQU5+8vDyNHDlSNWrUUEBAgHr37q2MjAynPsePH1dsbKwqV66s2rVrKz4+XgUFBU59Nm3apA4dOsjPz0/NmjXTggULrsrz1ltvqVGjRvL391dkZKS2b9/uchYAAACgInK5sOjevbsyMzMlSYGBgQoMDJQkZWdnq3v37i4da/PmzRo5cqS+/vprrV+/XpcuXVJ0dLRyc3MdfcaPH68VK1Zo6dKl2rx5s1JTU/XQQw852gsLCxUbG6uLFy9q69atWrhwoRYsWKDJkyc7+hw9elSxsbHq1q2bUlJSNG7cOD3xxBNau3ato8/ixYsVFxenKVOmaOfOnWrXrp1iYmJ06tSpYmcBAAAAKiqbMca48ga73a709PSrHtQ+deqU6tatq0uXLv3uMKdPn1bt2rW1efNm3XnnncrKylKtWrW0aNEi9enTR5J08OBBtWzZUklJSercubNWr16t++67T6mpqQoJCZEkzZs3TxMnTtTp06fl6+uriRMnatWqVdq7d6/jXP369VNmZqbWrFkjSYqMjFTHjh01Z84cSVJRUZHq16+v0aNH69lnny1WlhvJzs5WcHCwsrKyFBQU9LvHCcDVbNNcu2JaUZkpLn3kAwAqOFe+vxb7isXu3bu1e/duSdL+/fsd27t379a3336r9957T3Xr1nUreFZWliSpevXqkqTk5GRdunRJPXr0cPRp0aKFGjRooKSkJElSUlKS2rRp4ygqJCkmJkbZ2dnat2+fo8+Vx7jc5/IxLl68qOTkZKc+drtdPXr0cPQpThYAAACgoir2w9vh4eGy2Wyy2WzXvOWpUqVKevPNN393kKKiIo0bN0533HGHWrduLUlKT0+Xr6+vqlat6tQ3JCRE6enpjj5XFhWX2y+3/Vaf7Oxs/fzzzzp37pwKCwuv2efgwYPFzvJr+fn5ys/Pd2xnZ2ffaBgAAACAcqnYhcXRo0dljFGTJk20fft21apVy9Hm6+ur2rVry8vL63cHGTlypPbu3astW7b87mPcbKZPn65p06ZZHQMAAAAodcUuLBo2bCjplysLJW3UqFFauXKlvvjiC9WrV8+xPzQ0VBcvXlRmZqbTlYKMjAyFhoY6+vx69qbLMzVd2efXszdlZGQoKChIlSpVkpeXl7y8vK7Z58pj3CjLryUkJDit7ZGdna369esXZ0gAAACAcsXlWaEk6fDhwxo9erR69OihHj16aMyYMTp8+LDLxzHGaNSoUfrkk0+0YcMGNW7c2Kk9IiJCPj4+SkxMdOw7dOiQjh8/rqioKElSVFSU9uzZ4zR70/r16xUUFKRWrVo5+lx5jMt9Lh/D19dXERERTn2KioqUmJjo6FOcLL/m5+enoKAgpxcAAADgiVxeIG/t2rV64IEHFB4erjvuuEPSLytv33rrrVqxYoV69uxZ7GONHDlSixYt0meffabAwEDHswrBwcGqVKmSgoODNXToUMXFxal69eoKCgrS6NGjFRUV5ZiFKTo6Wq1atdKAAQM0c+ZMpaen6/nnn9fIkSPl5+cnSRo+fLjmzJmjCRMmaMiQIdqwYYOWLFmiVatWObLExcVp0KBBuu2229SpUyfNmjVLubm5evzxxx2ZbpQFAAAAqKhcnm62ffv2iomJ0YwZM5z2P/vss1q3bp127txZ/JNfZ0G9+fPna/DgwZJ+WZTu6aef1ocffqj8/HzFxMRo7ty5TrcfHTt2TCNGjNCmTZtUpUoVDRo0SDNmzJC39/+rmzZt2qTx48dr//79qlevniZNmuQ4x2Vz5szRK6+8ovT0dIWHh+uNN95QZGSko704WX4L080CpYfpZouH6WYBAK5w5fury4WFv7+/9uzZo1tuucVp/3fffae2bdsqLy/P9cQVBIUFUHooLIqHwgIA4IpSWcfislq1aiklJeWq/SkpKVctmgcAAACgYnD5GYthw4bpySef1JEjR3T77bdL+uUZi5dfftlpBiQAAAAAFYfLhcWkSZMUGBioV199VQkJCZKkOnXqaOrUqRozZkyJBwQAAABw83P5GYsrnT9/XpIUGBhYYoE8Gc9YAKWHZyyKh2csAACucOX7q8tXLK5EQQEAAABA+h0Pb2dkZGjAgAGqU6eOvL29HatWX34BAAAAqHhcvmIxePBgHT9+XJMmTVJYWNh116IAAAAAUHG4XFhs2bJFX375pcLDw0shDgAAAIDyyOVboerXry83nvcGAAAA4IFcLixmzZqlZ599Vj/++GMpxAEAAABQHhXrVqhq1ao5PUuRm5urpk2bqnLlyvLx8XHqe/bs2ZJNCAAAAOCmV6zCYtasWaUcAwAAAEB5VqzCYtCgQaWdAwAAAEA5VuxnLIqKivTyyy/rjjvuUMeOHfXss8/q559/Ls1sAAAAAMqJYhcWL730kp577jkFBASobt26mj17tkaOHFma2QAAAACUE8UuLD744APNnTtXa9eu1aeffqoVK1bo3//+t4qKikozHwAAAIByoNiFxfHjx3Xvvfc6tnv06CGbzabU1NRSCQYAAACg/Ch2YVFQUCB/f3+nfT4+Prp06VKJhwIAAABQvhRrVihJMsZo8ODB8vPzc+zLy8vT8OHDVaVKFce+jz/+uGQTAgAAALjpFbuwuNaUs4899liJhgEAAABQPhW7sJg/f35p5gAAAABQjhX7GQsAAAAAuB4KCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABus7Sw+OKLL3T//ferTp06stls+vTTT53ajTGaPHmywsLCVKlSJfXo0UPff/+9U5+zZ8/q0UcfVVBQkKpWraqhQ4cqJyfHqc/u3bv1xz/+Uf7+/qpfv75mzpx5VZalS5eqRYsW8vf3V5s2bfSf//zH5SwAAABARWVpYZGbm6t27drprbfeumb7zJkz9cYbb2jevHnatm2bqlSpopiYGOXl5Tn6PProo9q3b5/Wr1+vlStX6osvvtCTTz7paM/OzlZ0dLQaNmyo5ORkvfLKK5o6dar+8Y9/OPps3bpVjzzyiIYOHapvv/1WvXr1Uq9evbR3716XsgAAAAAVlc0YY6wOIUk2m02ffPKJevXqJemXKwR16tTR008/rWeeeUaSlJWVpZCQEC1YsED9+vXTgQMH1KpVK33zzTe67bbbJElr1qzRvffeq5MnT6pOnTp6++239de//lXp6eny9fWVJD377LP69NNPdfDgQUnSww8/rNzcXK1cudKRp3PnzgoPD9e8efOKlaU4srOzFRwcrKysLAUFBZXIuAH4hW2azeoI5YKZclN85AMAyglXvr/etM9YHD16VOnp6erRo4djX3BwsCIjI5WUlCRJSkpKUtWqVR1FhST16NFDdrtd27Ztc/S58847HUWFJMXExOjQoUM6d+6co8+V57nc5/J5ipPlWvLz85Wdne30AgAAADzRTVtYpKenS5JCQkKc9oeEhDja0tPTVbt2bad2b29vVa9e3anPtY5x5Tmu1+fK9htluZbp06crODjY8apfv/4N/q0BAACA8ummLSw8QUJCgrKyshyvEydOWB0JAAAAKBU3bWERGhoqScrIyHDan5GR4WgLDQ3VqVOnnNoLCgp09uxZpz7XOsaV57henyvbb5TlWvz8/BQUFOT0AgAAADzRTVtYNG7cWKGhoUpMTHTsy87O1rZt2xQVFSVJioqKUmZmppKTkx19NmzYoKKiIkVGRjr6fPHFF7p06ZKjz/r169W8eXNVq1bN0efK81zuc/k8xckCAAAAVGSWFhY5OTlKSUlRSkqKpF8ekk5JSdHx48dls9k0btw4vfjii1q+fLn27NmjgQMHqk6dOo6Zo1q2bKm7775bw4YN0/bt2/XVV19p1KhR6tevn+rUqSNJ6t+/v3x9fTV06FDt27dPixcv1uzZsxUXF+fIMXbsWK1Zs0avvvqqDh48qKlTp2rHjh0aNWqUJBUrCwAAAFCReVt58h07dqhbt26O7ctf9gcNGqQFCxZowoQJys3N1ZNPPqnMzEx16dJFa9askb+/v+M9//73vzVq1Cjdddddstvt6t27t9544w1He3BwsNatW6eRI0cqIiJCNWvW1OTJk53Wurj99tu1aNEiPf/883ruued0yy236NNPP1Xr1q0dfYqTBQAAAKiobpp1LCoC1rEASg/rWBQP61gAAFzhEetYAAAAACg/KCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuM3b6gAAAOC32abZrI5QLpgpxuoIQIXGFQsAAAAAbqOwAAAAAOA2CgsAAAAAbqOwAAAAAOA2Ht4GALiMh4mLh4eJAVQkXLEAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuo7AAAAAA4DYKCwAAAABuY7pZAACAKzCdcvEwnTJ+jSsWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANzGytsVCCuJFg8riQIAALiOwgIAAACW4YfP4ikPP3xyKxQAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhYWL3nrrLTVq1Ej+/v6KjIzU9u3brY4EAAAAWI5ZoVywePFixcXFad68eYqMjNSsWbMUExOjQ4cOqXbt2lbHw02GWS6KpzzMcgEAAG6MKxYueO211zRs2DA9/vjjatWqlebNm6fKlSvr/ffftzoaAAAAYCmuWBTTxYsXlZycrISEBMc+u92uHj16KCkp6Zrvyc/PV35+vmM7KytLkpSdnV26Ya8nz5rTljcl9t+H8S4WxrtsMd5li/EuW4x32WK8y5ZV3x8vn9eYG99hYDPF6QWlpqaqbt262rp1q6Kiohz7J0yYoM2bN2vbtm1XvWfq1KmaNm1aWcYEAAAAStyJEydUr1693+zDFYtSlJCQoLi4OMd2UVGRzp49qxo1ashm4/777Oxs1a9fXydOnFBQUJDVcTwe4122GO+yxXiXLca7bDHeZYvxdmaM0fnz51WnTp0b9qWwKKaaNWvKy8tLGRkZTvszMjIUGhp6zff4+fnJz8/PaV/VqlVLK2K5FRQUxP+4ZYjxLluMd9livMsW4122GO+yxXj/P8HBwcXqx8PbxeTr66uIiAglJiY69hUVFSkxMdHp1igAAACgIuKKhQvi4uI0aNAg3XbbberUqZNmzZql3NxcPf7441ZHAwAAACxFYeGChx9+WKdPn9bkyZOVnp6u8PBwrVmzRiEhIVZHK5f8/Pw0ZcqUq24XQ+lgvMsW4122GO+yxXiXLca7bDHevx+zQgEAAABwG89YAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt7GOBcrMmTNntHv3brVr107Vq1fXf//7X7333nvKz89X37591bJlS6sjAgAA4HdiHQuUie3btys6OlrZ2dmqWrWq1q9fr759+8rb21tFRUVKTU3Vli1b1KFDB6ujepwNGzZoy5YtSktLk91uV5MmTfTAAw/olltusTqaRyosLJSXl5dje9u2bcrPz1dUVJR8fHwsTOZZli1bpnvuuUeVK1e2OkqFdfToUf3www8KCwtT69atrY7jcXJzc5WcnOz02d2hQwfZbDaro3mkI0eOXPV3Zc+ePRUUFGR1tPLFAGWgR48e5oknnjDZ2dnmlVdeMfXq1TNPPPGEo/3xxx83vXr1sjCh58nIyDCdOnUydrvdeHt7G7vdbiIiIkxoaKjx8vIy8fHxVkf0KKmpqeaOO+4wXl5e5s477zRnz541sbGxxmazGZvNZv7whz+Y1NRUq2N6DJvNZoKCgsywYcPM119/bXUcjzdixAhz/vx5Y4wxFy5cML179zZ2u93YbDZjt9tNt27dHO1wT2FhoYmPjzeVK1c2drvdMc42m800bNjQLF++3OqIHiUnJ8f06dPHMcZ2u93x92RAQICZM2eO1RHLFZ6xQJlITk5WXFycAgMDNXbsWKWmpmrYsGGO9lGjRumbb76xMKHnGTNmjOrUqaNz584pJydHTz31lG699ValpaVp3bp1ev/99zV79myrY3qMiRMnyhijTz75RGFhYbrvvvuUnZ2tEydO6Mcff1StWrX00ksvWR3TozzzzDPasWOHoqKi1Lp1a82aNUtnzpyxOpZHeuedd3ThwgVJ0t/+9jdt27ZNn3/+uXJycvTFF1/o+PHj/PkuIc8995xWrlypxYsXa+3aterSpYtmzJih/fv3a+DAgerbt6/WrVtndUyPERcXp7S0NO3evVvfffedHnroIQ0cOFDZ2dmaPXu2JkyYoEWLFlkds/ywurJBxVClShVz9OhRx3ZAQIA5fPiwY/vYsWPG39/fgmSeKygoyOzdu9exnZOTY3x8fExWVpYxxpj/+Z//Mc2bN7cqnscJCwszSUlJxhhjzpw5Y2w2m/n8888d7YmJiaZJkyZWxfM4NpvNZGRkGGOM2bFjhxkxYoSpWrWq8fPzM3379jXr1q2zOKFnuXK8W7dubRYtWuTU/tlnn5k//OEPVkTzOGFhYeaLL75wbJ88edIEBASYvLw8Y4wxL7zwgomKirIqnsepWbOm2bFjh2P77Nmzxt/f3+Tm5hpjjJkzZ44JDw+3Kl65wxULlIn69evryJEjju2PPvpIYWFhju20tDTVrFnTimgey8/Pz+leXLvdrsLCQhUUFEiSbr/9dv34448WpfM8586dU926dSVJ1atXV+XKldWwYUNHe7NmzZSWlmZVPI8WERGhuXPnKi0tTe+++65Onz6tu+++W40bN7Y6mke5/HmSnp6utm3bOrW1a9dOJ06csCKWx8nJyXF8lkhSWFiY8vLydO7cOUlS7969tWvXLqvieZyCggKn5ygCAgJUUFCg3NxcSVJ0dLQOHjxoVbxyh8ICZaJfv346deqUYzs2NlaVKlVybC9fvlydOnWyIprH6tKliyZPnqzc3FxdunRJzz33nJo0aaLq1atLkk6fPq1q1apZnNJz1K5d26lwGDVqlGOspV8KjypVqlgRzSNd6wFWf39/DRgwQBs3btShQ4fUv39/C5J5rkmTJikuLk52u12pqalObWfOnOHPdwlp06aNPvzwQ8f2kiVLFBAQoNDQUElSUVGR/Pz8rIrncTp27Oh0W/Ds2bNVq1Yt1apVS9IvhV5AQIBV8codpptFmZgyZcpvtv/1r391mkkH7vv73/+u6OhoVa1aVTabTVWqVNHSpUsd7QcOHNDgwYOtC+hhwsPDlZSU5CiQZ8yY4dS+ZcuWq37lxe9nbjChYbNmzbjnvwTdeeedOnTokCSpVatWOnbsmFP7f/7zH916661WRPM4L7zwgmJjY7V8+XL5+/tr69ateuWVVxzta9asUfv27S1M6FlmzJihnj17atmyZfL19VV6eroWLlzoaN+6davuvfdeCxOWL0w3C3iwCxcu6KuvvlJ+fr46d+7M7WYW2r59uypXrsy0nCXk2LFjatCgAVNv3iSOHDkiX19f1atXz+ooHmHXrl1asmSJ8vPzFRMTo549e1odyaOlpaVp5cqVys/PV/fu3dWqVSurI5VbFBa4KXz22WfKysrSwIEDrY4CAACA34HCAjeFFi1a6Pvvv1dhYaHVUSqMHTt26MKFC7rzzjutjlIhpKWl6dKlS2rQoIHVUSoE/nyXLcYbnorPbtdQWAAVVMuWLfXdd99RzJURxrtsMd5li/EuO4x12WK8XcPD20AFlZiYqEuXLlkdo8L44IMPHAuMofTx57tsMd5lZ/r06crKyrI6RoXBZ7druGIBS2RmZmrp0qU6fvy4GjZsqL59+yo4ONjqWAAAAPidKCxQJh566CH1799fffr00b59+9S1a1fZbDY1adJEP/74o2w2mzZs2KCWLVtaHdXjpKena9u2bUpPT5ckhYaGKjIy0jEnOkrXpUuX5OPjY3UMj5Wbm6vk5GSlpaXJbrerSZMm6tChA7NFlbKCggJt3LjR8eNQt27dmDK8lC1YsEB//vOf+RGuFB0/ftzps6RGjRpWRyp/LFnvGxVOtWrVzIEDB4wxxtxzzz2mf//+Jj8/3xhjzMWLF83QoUNNdHS0lRE9Tk5Ojnn00UeNl5eX8fb2NrVr1za1a9c23t7exsvLyzz22GMmNzfX6pgeY/HixY4/08YY8+abb5oGDRoYu91uatSoYaZNm2ZhOs9TWFho4uPjTeXKlY3dbjd2u93YbDZjs9lMw4YNzfLly62O6FFGjRplVqxYYYwx5sSJE6ZFixbGy8vLhISEGC8vL9OmTRtz8uRJi1N6Nh8fH7N//36rY3ikt956y/F5feXrjjvuMDt27LA6XrnCytsoE3l5eY5fbVNSUvTMM8/I19dXkuTj46MJEyZo27ZtVkb0OGPHjtX27du1atUq5eXlKSMjQxkZGcrLy9N//vMfbd++XWPHjrU6psd45JFHlJmZKUmaP3++4uPjNXjwYK1YsULjx4/XzJkz9c9//tPakB7kueee08qVK7V48WKtXbtWXbp00YwZM7R//34NHDhQffv21bp166yO6TGWLl2qRo0aSZKefvpp1atXT+np6UpPT9epU6fUsGFDjRs3ztKMnqJ69erXfBUUFCgqKsqxjZLx97//XS+99JLi4+P1zjvvqHnz5po6dapWrVqlJk2a6M4779SOHTusjllucCsUykTnzp01dOhQDRs2TB06dNDkyZPVq1cvR/v69es1cOBApaWlWRfSw1SrVk2rVq3S7bfffs32r776Svfdd5/OnTtXxsk8k91uV3p6umrXrq3IyEj16dNH8fHxjva3335b7777rnbu3GlhSs9Rp04dLV68WH/84x8lST/99JNatGih//73v/Lz89Pf/vY3rV69Wlu3brU4qWeoVKmS9u/fr8aNG6t+/fpatmyZY5V5Sdq7d6+6deum06dPW5jSMwQGBupPf/qT+vbt69hnjNETTzyhF154QXXr1pUkDRo0yKqIHqVx48aaO3eu7rnnHknSd999p9tvv13p6eny9vbW2LFjdeDAAX6oKCZmhUKZmDRpkgYOHCgfHx+NGTNG48eP15kzZ9SyZUsdOnRIU6ZM0YABA6yO6VGKioocV4WuxdfXV0VFRWWYyPNdvq//yJEjio6OdmqLjo7WxIkTrYjlkXJychxfsCQpLCxMeXl5OnfunEJDQ9W7d2/NmDHDwoSe5Q9/+IO2b9+uxo0bKzAwUNnZ2U7t58+f5/OkhHz77bfq37+/NmzYoLfeeksBAQGSpGHDhqlXr16sCl3CTp065fR85y233KKsrCydPn1aYWFhGjJkiLp06WJhwvKFW6FQJmJjY/WPf/xDkyZN0tChQ3Xs2DENGzZMXbp00VNPPaXevXtr+vTpVsf0KPfdd5+efPJJffvtt1e1ffvttxoxYoTuv/9+C5J5rjVr1mj58uXy9/e/anrCvLw8HiguQW3atNGHH37o2F6yZIkCAgIckxIUFRXJz8/PqngeZ/z48XrmmWe0adMmJSQkaMyYMUpMTFRqaqo2btyov/zlL3rooYesjukRmjVrpq1btyo0NFTh4eH66quvrI7k0f7whz9o/fr1ju2NGzfK19fX8Vni7+/PZ7cLuGKBMtO7d2/16tVLycnJOnr0qIqKihQWFqaIiAgFBgZaHc/jzJkzR/3791dERISqVaum2rVrS/rl15nMzEzFxMRozpw5Fqf0LFfemrBhwwZFRUU5tr/++ms1bdrUilge6YUXXlBsbKyjkNu6dateeeUVR/uaNWvUvn17CxN6lsGDB+vs2bOKjY2VMUaFhYVOV+UeeOABvf766xYm9Cze3t56+eWXFRMTo/79++vRRx/ly20pSUhI0GOPPabPP/9c/v7++vjjjzVmzBjHeG/atEmtW7e2OGX5wTMWgIc7cOCAvv76a6fpZqOiotSiRQuLk1UsK1eulI+Pj2JiYqyO4jF27dqlJUuWKD8/XzExMerZs6fVkTxeZmam1q9fryNHjjh+HLrjjjt0yy23WB3NY505c0bDhg3Txo0b9fXXX6t58+ZWR/I4q1ev1r/+9S/HZ8mwYcMcbWfOnJEkpp4tJgoLlKkNGzZoy5YtTvNEP/DAA/ylBAAAUM5RWKBMnDp1Svfff7927Nghu92uoqIitW/fXj/99JNOnz6tuLg4zZw50+qYHoli7uZweSG3O++80+ooHuXIkSNX/fnu2bOngoKCrI7m0TIzM7V06VLHAnl9+/Zl4bZScuVihI0aNVLXrl1ZjBA3L6sW0EDF8vDDD5tevXqZrKwsk5eXZ0aNGmUGDhxojDEmMTHR1KhRw8yaNcvilJ4lIyPDdOrUydjtduPt7W3sdruJiIgwoaGhxsvLy8THx1sdsUJJSUkxdrvd6hgeIycnx/Tp08exKJ7dbnf82Q4ICDBz5syxOqJH+fOf/2yWLl1qjDFm7969pmbNmqZWrVomMjLShISEmNDQUBZvKyEsRli2Ll68aOLj403Tpk1Nx44dzXvvvefUnp6ezme3C5gVCmVi9erVevHFFxUUFCQ/Pz/NmDFDH374obKzs9W9e3fNmjVLb7/9ttUxPcqYMWNUp04dnTt3Tjk5OXrqqad06623Ki0tTevWrdP777+v2bNnWx0T+F3i4uKUlpam3bt367vvvtNDDz2kgQMHKjs7W7Nnz9aECRO0aNEiq2N6jCsfYI2Pj1d0dLROnjypr7/+WidOnFBsbCwL5JUQFiMsWy+99JI++OADDR8+XNHR0YqLi9Nf/vIXpz6Gm3uKz+rKBhVDrVq1zL59+xzbFy5cMHa73Zw5c8YYY8zhw4eNn5+fVfE8UlBQkNm7d69jOycnx/j4+JisrCxjjDH/8z//Y5o3b25VPI9TrVq133wFBQXxq1cJqlmzptmxY4dj++zZs8bf39/k5uYaY4yZM2eOCQ8Ptyqex6lUqZL54YcfjDHGhIWFmZ07dzq1Hzp0yAQHB1uQzPP4+/ubI0eOGGOMqVevntm2bZtT+549e0zNmjWtiOaRmjVr5rhCZIwx33//vWnWrJkZPHiwKSoq4oqFi5huFmWiS5cumjx5shYuXChfX18999xzatKkiapXry5JOn36tKpVq2ZxSs/i5+fnND2h3W5XYWGhCgoKJEm33367fvzxR4vSeZ78/HyNGDFCbdq0uWb7sWPHNG3atDJO5bkKCgqcnqMICAhQQUGBcnNzVblyZUVHR+uZZ56xMKFnadu2rTZs2KCmTZsqNDRUx44dc5rO99ixY6pUqZKFCT0HixGWrZ9++slpOtlmzZpp06ZN6t69uwYMGMDzny6isECZ+Pvf/67o6GhVrVpVNptNVapU0dKlSx3tBw4c0ODBg60L6IEo5spWeHi46tev77SWxZV27dpFYVGCOnbsqNmzZzvWYpk9e7Zq1aqlWrVqSfplZe7LKxbDfZMmTdLAgQPl4+OjMWPGaPz48Tpz5oxatmypQ4cOacqUKRowYIDVMT3C5cUIQ0JCHIsRvvnmm46xHjt2LIsRlqDQ0FAdPnzYcfuZJNWtW1cbN25Ut27d+G7iIgoLlIkmTZpo9+7d2rJliy5evKjOnTurZs2aMsbIZrPxP24poJgrW7GxscrMzLxue/Xq1TVw4MCyC+ThZsyYoZ49e2rZsmXy9fVVenq6Fi5c6GjfunWr7r33XgsTepbY2Fj94x//0Lhx45SamipjjGOufz8/Pw0fPlzTp0+3OKVnYDHCstW9e3ctWrRId911l9P+OnXqaMOGDeratas1wcopppuFpXx9fbVr1y61bNnS6ige6cKFC/rqq6+Un5/vKOYAT5GWlqaVK1cqPz9f3bt3V6tWrayO5PEKCwu1c+dOpwXyIiIiFBgYaHU0j5OZmal169bp6NGjLEZYio4dO6aDBw9ed/HS1NRUrV+//rpXo+GMwgJlIi4u7pr7Z8+erccee8yxouVrr71WlrEAAABQQrgVCmVi1qxZateunapWreq03xijAwcOqEqVKk4PGqP0ZWRk6J133tHkyZOtjuJRTp48qapVq151f/+lS5eUlJTEAnkl6MyZM9q9e7fatWun6tWr67///a/ee+895efnq2/fvlwJLUXGGG3atEk//PCDwsLCFBMTIx8fH6tjeZTt27crKSlJ6enpkn55FiAqKkqdOnWyOJlnyc/Pl91ud/z5PXz4sN5//33H4o9Dhw5V48aNLU5Zjlg2HxUqlOnTp5vGjRubxMREp/3e3t5O09Ci7LBgW8lKTU01HTt2NHa73Xh5eZkBAwaY8+fPO9qZsrBkbdu2zQQHBxubzWaqVatmduzYYRo3bmxuueUW07RpU1OpUiWTnJxsdUyPcc8995jMzExjjDFnzpwxkZGRxmazmVq1ahm73W5atGhhTp06ZXFKz5CRkWG6dOlibDabadiwoenUqZPp1KmTadiwobHZbKZLly4mIyPD6pge409/+pNj8cctW7YYPz8/07ZtW/Pwww+b9u3bm8qVK5utW7danLL84FYolJlvvvlGjz32mO6//35Nnz5dPj4+8vHx0a5du7g3uhTs3r37N9sPHjyoRx55RIWFhWWUyLMNGjRIhw4d0pw5c5SZmalnn31WNptN69atU7Vq1ZSRkaGwsDCmiSwhPXv2VKNGjfTaa6/pnXfe0ezZs3X33Xfr3XfflSQNGTJE586d0yeffGJxUs9gt9uVnp6u2rVr66mnntLmzZu1cuVKNW7cWCdPnlSvXr3UsWNHFjotAX369FFqaqrmz5+v5s2bO7UdOnRIQ4YMUZ06dZwm48DvFxwcrB07duiWW25R165d1aFDB6fbsidNmqSNGzdqy5YtFqYsR6yubFCxnD9/3gwcONC0bdvW7Nmzx/j4+HDFopTYbDZjt9uNzWa76nV5P7+gl5w6deo4LWSVl5dn7r//fhMeHm7OnDnDFYsSVq1aNbN//35jjDEXL140drvdafyTk5NN3bp1rYrncWw2m+NX8ubNm5vPPvvMqf3zzz83jRs3tiKaxwkICLhqAcIr7dixwwQEBJRhIs9WpUoVc+DAAWOMMSEhISYlJcWp/YcffmC8XWC3urBBxRIQEKCFCxcqISFBPXr04NfyUlS9enW9++67Onr06FWvI0eOaOXKlVZH9ChZWVlO64L4+fnp448/VqNGjdStWzedOnXKwnSe5+LFi44F2Xx8fFS5cmWnWc9q1qypM2fOWBXPI11+Du7cuXNq2rSpU1uzZs2UmppqRSyP4+fnd9WieFc6f/68/Pz8yjCRZ4uMjNSKFSskSU2bNtWuXbuc2lNSUhzrP+HGeHgblujXr5+6dOmi5ORkNWzY0Oo4HikiIkKpqanXHd/MzEwZ7oQsMZfXarlyKkhvb28tXbpUffv21X333WdhOs9Tv359HTlyxLGo1UcffaSwsDBHe1paGtMrl7DBgwfLz89Ply5d0tGjR3Xrrbc62tLT06+anAO/z8MPP6xBgwbp9ddf11133eVYYT47O1uJiYmKi4vTI488YnFKz/Hiiy/qnnvuUW5urh555BE9/fTT+v777x0LEr7xxhtKSEiwOma5QWEBy9SrV0/16tWzOobHGj58uHJzc6/b3qBBA82fP78ME3m2e+65R//4xz/Uu3dvp/2Xi4vevXvrxIkTFqXzPP369XO6ChQbG+vUvnz5cmbPKUFXzuH/4IMP6sKFC07ty5YtU3h4eBmn8kyvvfaaioqK1K9fPxUUFMjX11fSL1fpvL29NXToUP3973+3OKXniIqK0urVqxUXF6dt27ZJkl566SVJvyySN3XqVI0dO9bKiOUKD28DQAkoKCjQhQsXHL8uXqv9p59+4gpdGblw4YK8vLy4ZaSM5ObmysvLS/7+/lZH8RjZ2dlKTk52mm42IiLiup8xcN/p06edFn+8fEUUxcczFkAFdeLECQ0ZMsTqGB7D29v7N//CT0tL07Rp08owUcV25swZjRgxwuoYFcbZs2f11FNPWR3DYxw4cEDLli1TWFiYHnnkEbVv315LlizRuHHjtGHDBqvjeZwDBw5o/vz5Onv2rCIjI1WtWjW9/PLLGjJkCOPtIq5YABXUrl271KFDBx6gLyOMd9livMsW411y1qxZowcffFABAQG6cOGCPvnkEw0cOFDt2rVTUVGRNm/erHXr1ql79+5WR/UIjHfJ4hkLwEMtX778N9uPHDlSRkkqBsa7bDHeZYvxLjsvvPCC4uPj9eKLL+qjjz5S//79NWLECMd9/wkJCZoxYwZfdEsI412yuGIBeCi73S6bzfabMz/ZbDZ+YSwhjHfZYrzLFuNddoKDg5WcnKxmzZqpqKhIfn5+2r59u9q3by9J2rt3r3r06OF49gLuYbxLFs9YAB4qLCxMH3/8sYqKiq752rlzp9URPQrjXbYY77LFeJety2uG2O12+fv7Kzg42NEWGBiorKwsq6J5JMa75FBYAB4qIiJCycnJ122/0a+PcA3jXbYY77LFeJedRo0a6fvvv3dsJyUlqUGDBo7t48ePO63ZAvcw3iWLZywADxUfH/+b61g0a9ZMGzduLMNEno3xLluMd9livMvOiBEjnG4pa926tVP76tWrud+/BDHeJYtnLAAAAAC4jVuhAAAAALiNwgIAAACA2ygsAAAAALiNwgIAAACA2ygsAAAAALiNwgIAgJvMn//8Z1WrVk19+vSxOgoAFBuFBQAAN5mxY8fqgw8+sDoGALiEwgIAgP/fmTNnVLt2bf3444+W5ujatasCAwOv2t+vXz+9+uqrFiQCgBujsAAADzR48GDZbLarXj/88IPV0W5qL730kh588EE1atTI6ijX9Pzzz+ull15SVlaW1VEA4CreVgcAAJSOu+++W/Pnz3faV6tWrav6Xbx4Ub6+vmUV66Z14cIFvffee1q7dm2pnys8PFwFBQVX7V+3bp3q1Klz3fe1bt1aTZs21b/+9S+NHDmyNCMCgMu4YgEAHsrPz0+hoaFOLy8vL3Xt2lWjRo3SuHHjVLNmTcXExEiSioqKNH36dDVu3FiVKlVSu3bt9L//+79Ox8zNzdXAgQMVEBCgsLAwvfrqq+ratavGjRvn6NOoUSPNmjXL6X3h4eGaOnVqsc/TtWtXjRkzRhMmTFD16tUVGhrqeP9lRUVFmjlzppo1ayY/Pz81aNBAL730kj744APVqFFD+fn5Tv179eqlAQMGXHe8/vOf/8jPz0+dO3e+4TmuzDl69GiNGzdO1apVU0hIiN59913l5ubq8ccfV2BgoJo1a6bVq1c7nSslJUV79+696vVbRcVl999/vz766KMb9gOAskZhAQAV0MKFC+Xr66uvvvpK8+bNkyRNnz5dH3zwgebNm6d9+/Zp/Pjxeuyxx7R582bH++Lj47V582Z99tlnWrdunTZt2qSdO3e6dO7inOdyxipVqmjbtm2aOXOmXnjhBa1fv97RnpCQoBkzZmjSpEnav3+/Fi1apJCQEPXt21eFhYVavny5o++pU6e0atUqDRky5Lq5vvzyS0VERDjtu945fp2zZs2a2r59u0aPHq0RI0aob9++uv3227Vz505FR0drwIABunDhgkvjdD2dOnXS9u3bryqcAMByBgDgcQYNGmS8vLxMlSpVHK8+ffoYY4z505/+ZNq3b+/UPy8vz1SuXNls3brVaf/QoUPNI488Yowx5vz588bX19csWbLE0X7mzBlTqVIlM3bsWMe+hg0bmtdff93pOO3atTNTpkwp1nkuZ+zSpYtTn44dO5qJEycaY4zJzs42fn5+5t13373mv/+IESPMPffc49h+9dVXTZMmTUxRUdE1+xtjzIMPPmiGDBni2L7ROa6Vs6CgwFSpUsUMGDDAsS8tLc1IMklJSdc9zq/dddddpmbNmqZSpUqmbt26TuO1a9cuI8n8+OOPxT4eAJQFnrEAAA/VrVs3vf32247tKlWqOP7517/M//DDD7pw4YJ69uzptP/ixYtq3769JOnw4cO6ePGiIiMjHe3Vq1dX8+bNi52pOOe5rG3btk7bYWFhOnXqlCTpwIEDys/P11133XXN8wwbNkwdO3bUTz/9pLp162rBggWOB9qv5+eff5a/v79j+0bnuFZOLy8v1ahRQ23atHHsu3yF43L24vj888+v21apUiVJKrErIABQUigsAMBDValSRc2aNbtu25VycnIkSatWrVLdunWd2vz8/Fw6r91ulzHGad+lS5dcPo+Pj4/Tts1mU1FRkaT/9+X6etq3b6927drpgw8+UHR0tPbt26dVq1b95ntq1qypc+fOObZvdI7fynnlvsvFzOXs7jp79qykaz+IDwBW4hkLAIBatWolPz8/HT9+XM2aNXN61a9fX5LUtGlT+fj4aNu2bY73nTt3Tt99953TsWrVqqW0tDTHdnZ2to4ePVrs8xTHLbfcokqVKikxMfG6fZ544gktWLBA8+fPV48ePW54/Pbt22v//v0uncMKe/fuVb169VSzZk2rowCAE65YAAAUGBioZ555RuPHj1dRUZG6dOmirKwsffXVVwoKCtKgQYMUEBCgoUOHKj4+XjVq1FDt2rX117/+VXa7829U3bt314IFC3T//feratWqmjx5sry8vIp9nuLw9/fXxIkTNWHCBPn6+uqOO+7Q6dOntW/fPg0dOlSS1L9/fz3zzDN69913i7WKdUxMjBISEnTu3DlVq1atWOewwpdffqno6GjLzg8A10NhAQCQJP3tb39TrVq1NH36dB05ckRVq1ZVhw4d9Nxzzzn6vPLKK8rJydH999+vwMBAPf3001ct1paQkKCjR4/qvvvuU3BwsP72t785rlgU9zzFMWnSJHl7e2vy5MlKTU1VWFiYhg8f7mgPDg5W7969tWrVKvXq1euGx2vTpo06dOigJUuW6C9/+UuxzlHW8vLy9Omnn2rNmjWWZQCA67GZX98ICwCAC7p27arw8PCr1q64Gdx111269dZb9cYbbxSr/6pVqxQfH6+9e/dedSXmZvD222/rk08+0bp166yOAgBX4YoFAMDjnDt3Tps2bdKmTZs0d+7cYr8vNjZW33//vX766SeXnvkoKz4+PnrzzTetjgEA10RhAQDwOO3bt9e5c+f08ssvuzQdriSnVcRvNk888YTVEQDgurgVCgAAAIDbbr4bSAEAAACUOxQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbf8f8cqvi++vc+EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample2 = []\n",
    "for ele in sample.items():\n",
    "    # print(ele[0].state)\n",
    "    sample2.append(ele[0].state * ele[1])\n",
    "counts2 = np.sum(sample2, axis=0)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.ylabel(\"Photon number\")\n",
    "plt.xlabel(r\"Frequency (cm$^{-1}$)\")\n",
    "plt.xticks(range(len(omegap)), np.round(omegap, 1), rotation=90)\n",
    "plt.bar(range(len(omegap)), counts2, color=\"green\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在模拟一个光激发过程，其中涉及从电子基态的预激发振动态进行振动跃迁。\n",
    "预激发振动态可以通过应用位移门来模拟。\n",
    "我们向第6个振动模式插入平均一个振动量子，并计算激发电子态中每个振动模式的平均光子数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"6.363636363636363cm\" version=\"1.1\" width=\"14.7cm\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><polyline fill=\"none\" points=\"40,180 130,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"175\" /><text font-size=\"9\" x=\"83\" y=\"170\">D</text><text font-size=\"7\" x=\"95\" y=\"168\">r =1.0</text><text font-size=\"7\" x=\"95\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"40,30 130,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"25\" /><text font-size=\"9\" x=\"83\" y=\"20\">D</text><text font-size=\"7\" x=\"95\" y=\"18\">r =0.056</text><text font-size=\"7\" x=\"95\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"40,60 130,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"55\" /><text font-size=\"9\" x=\"83\" y=\"50\">D</text><text font-size=\"7\" x=\"95\" y=\"48\">r =-0.053</text><text font-size=\"7\" x=\"95\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"40,90 130,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"85\" /><text font-size=\"9\" x=\"83\" y=\"80\">D</text><text font-size=\"7\" x=\"95\" y=\"78\">r =0.982</text><text font-size=\"7\" x=\"95\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"40,120 130,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"115\" /><text font-size=\"9\" x=\"83\" y=\"110\">D</text><text font-size=\"7\" x=\"95\" y=\"108\">r =0.525</text><text font-size=\"7\" x=\"95\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"40,150 130,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"145\" /><text font-size=\"9\" x=\"83\" y=\"140\">D</text><text font-size=\"7\" x=\"95\" y=\"138\">r =-0.112</text><text font-size=\"7\" x=\"95\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"130,180 220,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"172.5\" y=\"175\" /><text font-size=\"9\" x=\"173\" y=\"170\">D</text><text font-size=\"7\" x=\"185\" y=\"168\">r =0.525</text><text font-size=\"7\" x=\"185\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"40,210 130,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"205\" /><text font-size=\"9\" x=\"83\" y=\"200\">D</text><text font-size=\"7\" x=\"95\" y=\"198\">r =-0.016</text><text font-size=\"7\" x=\"95\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"220,30 240,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,30 310,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,60 240,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,60 310,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,90 240,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,90 310,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,120 240,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,120 310,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,150 240,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,150 310,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,180 240,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,180 310,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"220,210 240,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"290,210 310,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"240\" y=\"20\" /><text font-size=\"10\" x=\"261\" y=\"120.0\">U</text><polyline fill=\"none\" points=\"310,30 400,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"25\" /><text font-size=\"9\" x=\"353\" y=\"20\">S</text><text font-size=\"7\" x=\"365\" y=\"18\">r =-0.097</text><text font-size=\"7\" x=\"365\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"310,60 400,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"55\" /><text font-size=\"9\" x=\"353\" y=\"50\">S</text><text font-size=\"7\" x=\"365\" y=\"48\">r =-0.07</text><text font-size=\"7\" x=\"365\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"310,90 400,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"85\" /><text font-size=\"9\" x=\"353\" y=\"80\">S</text><text font-size=\"7\" x=\"365\" y=\"78\">r =-0.021</text><text font-size=\"7\" x=\"365\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"310,120 400,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"115\" /><text font-size=\"9\" x=\"353\" y=\"110\">S</text><text font-size=\"7\" x=\"365\" y=\"108\">r =0.06</text><text font-size=\"7\" x=\"365\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"310,150 400,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"145\" /><text font-size=\"9\" x=\"353\" y=\"140\">S</text><text font-size=\"7\" x=\"365\" y=\"138\">r =0.075</text><text font-size=\"7\" x=\"365\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"310,180 400,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"175\" /><text font-size=\"9\" x=\"353\" y=\"170\">S</text><text font-size=\"7\" x=\"365\" y=\"168\">r =0.112</text><text font-size=\"7\" x=\"365\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"310,210 400,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"352.5\" y=\"205\" /><text font-size=\"9\" x=\"353\" y=\"200\">S</text><text font-size=\"7\" x=\"365\" y=\"198\">r =0.187</text><text font-size=\"7\" x=\"365\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"400,30 420,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,30 490,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,60 420,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,60 490,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,90 420,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,90 490,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,120 420,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,120 490,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,150 420,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,150 490,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,180 420,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,180 490,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"400,210 420,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"470,210 490,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"420\" y=\"20\" /><text font-size=\"10\" x=\"441\" y=\"120.0\">U</text><polyline fill=\"none\" points=\"130,30 220,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,60 220,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,90 220,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,120 220,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,150 220,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,210 220,210\" stroke=\"black\" stroke-width=\"2\" /><text font-size=\"12\" x=\"25\" y=\"30\">0</text><text font-size=\"12\" x=\"25\" y=\"60\">1</text><text font-size=\"12\" x=\"25\" y=\"90\">2</text><text font-size=\"12\" x=\"25\" y=\"120\">3</text><text font-size=\"12\" x=\"25\" y=\"150\">4</text><text font-size=\"12\" x=\"25\" y=\"180\">5</text><text font-size=\"12\" x=\"25\" y=\"210\">6</text></svg>"
      ],
      "text/plain": [
       "<svgwrite.drawing.Drawing at 0x155547680>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cir = dq.photonic.QumodeCircuit(\n",
    "    nmode=modes,\n",
    "    init_state=\"vac\",\n",
    "    # init_state=init_state,\n",
    "    cutoff=cutoff,\n",
    "    backend=\"gaussian\",\n",
    ")\n",
    "\n",
    "cir.d(wires=[5], r=1.0)\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.d(wires=[i], r=delta_2[i])\n",
    "\n",
    "cir.any(cr, wires=list(range(modes)))\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.s(wires=[i], r=-lambda_2[i])\n",
    "\n",
    "cir.any(cl, wires=list(range(modes)))\n",
    "\n",
    "state = cir()\n",
    "\n",
    "# 线路可视化\n",
    "cir.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "chain 1: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:03<00:00, 26099.72it/s]\u001b[0m\n",
      "chain 2: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 45578.62it/s]\u001b[0m\n",
      "chain 3: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 48048.44it/s]\u001b[0m\n",
      "chain 4: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 49628.73it/s]\u001b[0m\n",
      "chain 5: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 48353.82it/s]\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sample3 = cir.measure(shots=shots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGGCAYAAADmRxfNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQ3ElEQVR4nO3deXhU5f3//9dMyAJkYc0GYVcgsgQihCi1oJCocaEGLhBkkYiFsibKEquAVAtirWBRUdsKtqUCH5cikQAGQYEIEgRldQEFTSYgkIQEk5Dk/P7wx/kygpJxkhwyeT6ua656zn3POa/eTcd5zznnvm2GYRgCAAAAADfYrQ4AAAAAoPajsAAAAADgNgoLAAAAAG6jsAAAAADgNgoLAAAAAG6jsAAAAADgNgoLAAAAAG6jsAAAAADgtnpWB6hLKioqlJ2drYCAANlsNqvjAAAAAL/IMAydPXtW4eHhstt/+ZoEhUUNys7OVkREhNUxAAAAAJccP35cLVu2/MU+FBY1KCAgQNKP/8MEBgZanAYAAAD4ZQUFBYqIiDC/x/4SCosadOH2p8DAQAoLAAAA1BqVuY2fh7cBAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuI3CAgAAAIDbKCwAAAAAuK2e1QEAoCrYHrdZHaFWMOYYVkcAAHgorlgAAAAAcBuFBQAAAAC3UVgAAAAAcJvlhcV3332n++67T02bNlX9+vXVtWtX7dq1y2w3DEOzZ89WWFiY6tevrwEDBuiLL75wOsbp06c1YsQIBQYGqlGjRkpKSlJhYaFTn08//VS/+c1v5Ofnp4iICC1cuPCSLKtXr1anTp3k5+enrl276t1333Vqr0wWAAAAoC6ytLA4c+aMbrzxRnl7e2vdunU6cOCAnnnmGTVu3Njss3DhQj333HNaunSpduzYoYYNGyo+Pl7FxcVmnxEjRmj//v3auHGj1q5dqw8++EAPPvig2V5QUKC4uDi1bt1aWVlZevrppzV37ly9/PLLZp/t27fr3nvvVVJSkj755BMNGjRIgwYN0r59+1zKAgAAANRFNsMwLJsiZNasWdq2bZs+/PDDy7YbhqHw8HA99NBDevjhhyVJ+fn5CgkJ0bJlyzRs2DAdPHhQkZGR+vjjj3X99ddLktLT03X77bfr22+/VXh4uF588UX98Y9/lMPhkI+Pj3nut99+W4cOHZIkDR06VEVFRVq7dq15/j59+igqKkpLly6tVJYrKSgoUFBQkPLz8xUYGPjrBw7AJZgVqnKYFQoA4ApXvr9aesVizZo1uv766zVkyBAFBwerR48eeuWVV8z2o0ePyuFwaMCAAea+oKAgxcTEKDMzU5KUmZmpRo0amUWFJA0YMEB2u107duww+9x0001mUSFJ8fHxOnz4sM6cOWP2ufg8F/pcOE9lsvxUSUmJCgoKnF4AAACAJ7K0sDhy5IhefPFFXXPNNVq/fr0mTJigKVOmaPny5ZIkh8MhSQoJCXF6X0hIiNnmcDgUHBzs1F6vXj01adLEqc/ljnHxOX6uz8XtV8ryU/Pnz1dQUJD5ioiIuNKQAAAAALWSpYVFRUWFevbsqT//+c/q0aOHHnzwQY0bN05Lly61MlaVSU1NVX5+vvk6fvy41ZEAAACAamFpYREWFqbIyEinfZ07d9axY8ckSaGhoZKk3Nxcpz65ublmW2hoqE6cOOHUXlZWptOnTzv1udwxLj7Hz/W5uP1KWX7K19dXgYGBTi8AAADAE1laWNx44406fPiw077PP/9crVu3liS1bdtWoaGhysjIMNsLCgq0Y8cOxcbGSpJiY2OVl5enrKwss8+mTZtUUVGhmJgYs88HH3yg8+fPm302btyojh07mjNQxcbGOp3nQp8L56lMFgAAAKCusrSwSE5O1kcffaQ///nP+vLLL7VixQq9/PLLmjhxoiTJZrNp2rRpeuKJJ7RmzRp99tlnGjVqlMLDwzVo0CBJP17huPXWWzVu3Djt3LlT27Zt06RJkzRs2DCFh4dLkoYPHy4fHx8lJSVp//79WrlypRYvXqyUlBQzy9SpU5Wenq5nnnlGhw4d0ty5c7Vr1y5NmjSp0lkAAACAuqqelSfv1auX3nrrLaWmpmrevHlq27atFi1apBEjRph9ZsyYoaKiIj344IPKy8tT3759lZ6eLj8/P7PPf/7zH02aNEm33HKL7Ha7EhMT9dxzz5ntQUFB2rBhgyZOnKjo6Gg1a9ZMs2fPdlrr4oYbbtCKFSv06KOP6pFHHtE111yjt99+W126dHEpCwAAAFAXWbqORV3DOhZA9WEdi8phHQsAgCtqzToWAAAAADwDhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHCbpYXF3LlzZbPZnF6dOnUy24uLizVx4kQ1bdpU/v7+SkxMVG5urtMxjh07poSEBDVo0EDBwcGaPn26ysrKnPps3rxZPXv2lK+vrzp06KBly5ZdkuX5559XmzZt5Ofnp5iYGO3cudOpvTJZAAAAgLrK8isW1113nXJycszX1q1bzbbk5GS98847Wr16tbZs2aLs7Gzdc889Znt5ebkSEhJUWlqq7du3a/ny5Vq2bJlmz55t9jl69KgSEhLUv39/7dmzR9OmTdMDDzyg9evXm31WrlyplJQUzZkzR7t371b37t0VHx+vEydOVDoLAAAAUJfZDMMwrDr53Llz9fbbb2vPnj2XtOXn56t58+ZasWKFBg8eLEk6dOiQOnfurMzMTPXp00fr1q3THXfcoezsbIWEhEiSli5dqpkzZ+rkyZPy8fHRzJkzlZaWpn379pnHHjZsmPLy8pSeni5JiomJUa9evbRkyRJJUkVFhSIiIjR58mTNmjWrUlkqo6CgQEFBQcrPz1dgYOCvHjcAl7I9brM6Qq1gzLHsIx8AUAu58v3V8isWX3zxhcLDw9WuXTuNGDFCx44dkyRlZWXp/PnzGjBggNm3U6dOatWqlTIzMyVJmZmZ6tq1q1lUSFJ8fLwKCgq0f/9+s8/Fx7jQ58IxSktLlZWV5dTHbrdrwIABZp/KZLmckpISFRQUOL0AAAAAT2RpYRETE6Nly5YpPT1dL774oo4eParf/OY3Onv2rBwOh3x8fNSoUSOn94SEhMjhcEiSHA6HU1Fxof1C2y/1KSgo0A8//KDvv/9e5eXll+1z8TGulOVy5s+fr6CgIPMVERFRuYEBAAAAapl6Vp78tttuM/+5W7duiomJUevWrbVq1SrVr1/fwmRVIzU1VSkpKeZ2QUEBxQUAAAA8kuW3Ql2sUaNGuvbaa/Xll18qNDRUpaWlysvLc+qTm5ur0NBQSVJoaOglMzNd2L5Sn8DAQNWvX1/NmjWTl5fXZftcfIwrZbkcX19fBQYGOr0AAAAAT3RVFRaFhYX66quvFBYWpujoaHl7eysjI8NsP3z4sI4dO6bY2FhJUmxsrD777DOn2Zs2btyowMBARUZGmn0uPsaFPheO4ePjo+joaKc+FRUVysjIMPtUJgsAAABQl1l6K9TDDz+sO++8U61bt1Z2drbmzJkjLy8v3XvvvQoKClJSUpJSUlLUpEkTBQYGavLkyYqNjTVnYYqLi1NkZKRGjhyphQsXyuFw6NFHH9XEiRPl6+srSRo/fryWLFmiGTNmaOzYsdq0aZNWrVqltLQ0M0dKSopGjx6t66+/Xr1799aiRYtUVFSk+++/X5IqlQUAAACoyywtLL799lvde++9OnXqlJo3b66+ffvqo48+UvPmzSVJzz77rOx2uxITE1VSUqL4+Hi98MIL5vu9vLy0du1aTZgwQbGxsWrYsKFGjx6tefPmmX3atm2rtLQ0JScna/HixWrZsqX+/ve/Kz4+3uwzdOhQnTx5UrNnz5bD4VBUVJTS09OdHui+UhYAAACgLrN0HYu6hnUsgOrDOhaVwzoWAABX1Kp1LAAAAADUfhQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbS4VFufPn1f79u118ODB6soDAAAAoBZyqbDw9vZWcXFxdWUBAAAAUEu5fCvUxIkT9dRTT6msrKw68gAAAACoheq5+oaPP/5YGRkZ2rBhg7p27aqGDRs6tb/55ptVFg4AAABA7eByYdGoUSMlJiZWRxYAAAAAtZTLhcWrr75aHTkAAAAA1GK/arrZsrIyvffee3rppZd09uxZSVJ2drYKCwurNBwAAACA2sHlKxbffPONbr31Vh07dkwlJSUaOHCgAgIC9NRTT6mkpERLly6tjpwAAAAArmIuX7GYOnWqrr/+ep05c0b169c39//ud79TRkZGlYYDAAAAUDu4fMXiww8/1Pbt2+Xj4+O0v02bNvruu++qLBgAAACA2sPlKxYVFRUqLy+/ZP+3336rgICAKgkFAAAAoHZxubCIi4vTokWLzG2bzabCwkLNmTNHt99+e1VmAwAAAFBLuHwr1DPPPKP4+HhFRkaquLhYw4cP1xdffKFmzZrpv//9b3VkBAAAAHCVc7mwaNmypfbu3avXX39dn376qQoLC5WUlKQRI0Y4PcwNAAAAoO5wubCQpHr16um+++6r6iwAAAAAaqlftUDe4cOHNWnSJN1yyy265ZZbNGnSJB06dMitIAsWLJDNZtO0adPMfcXFxZo4caKaNm0qf39/JSYmKjc31+l9x44dU0JCgho0aKDg4GBNnz5dZWVlTn02b96snj17ytfXVx06dNCyZcsuOf/zzz+vNm3ayM/PTzExMdq5c6dTe2WyAAAAAHWVy4XFG2+8oS5duigrK0vdu3dX9+7dtXv3bnXt2lVvvPHGrwrx8ccf66WXXlK3bt2c9icnJ+udd97R6tWrtWXLFmVnZ+uee+4x28vLy5WQkKDS0lJt375dy5cv17JlyzR79myzz9GjR5WQkKD+/ftrz549mjZtmh544AGtX7/e7LNy5UqlpKRozpw52r17t7p37674+HidOHGi0lkAAACAusxmGIbhyhvat2+vESNGaN68eU7758yZo3//+9/66quvXApQWFionj176oUXXtATTzyhqKgoLVq0SPn5+WrevLlWrFihwYMHS5IOHTqkzp07KzMzU3369NG6det0xx13KDs7WyEhIZKkpUuXaubMmTp58qR8fHw0c+ZMpaWlad++feY5hw0bpry8PKWnp0uSYmJi1KtXLy1ZskTSj1PqRkREaPLkyZo1a1alslRGQUGBgoKClJ+fr8DAQJfGCcAvsz1uszpCrWDMcekjHwBQx7ny/dXlKxY5OTkaNWrUJfvvu+8+5eTkuHo4TZw4UQkJCRowYIDT/qysLJ0/f95pf6dOndSqVStlZmZKkjIzM9W1a1ezqJCk+Ph4FRQUaP/+/Wafnx47Pj7ePEZpaamysrKc+tjtdg0YMMDsU5ksAAAAQF3m8sPb/fr104cffqgOHTo47d+6dat+85vfuHSs119/Xbt379bHH398SZvD4ZCPj48aNWrktD8kJEQOh8Psc3FRcaH9Qtsv9SkoKNAPP/ygM2fOqLy8/LJ9Ljw3Upksl1NSUqKSkhJzu6Cg4Gf7AgAAALVZpQqLNWvWmP981113aebMmcrKyjJvAfroo4+0evVqPf7445U+8fHjxzV16lRt3LhRfn5+LsauHebPn+/SmAAAAAC1VaUKi0GDBl2y74UXXtALL7zgtG/ixIkaP358pU6clZWlEydOqGfPnua+8vJyffDBB1qyZInWr1+v0tJS5eXlOV0pyM3NVWhoqCQpNDT0ktmbLszUdHGfn87elJubq8DAQNWvX19eXl7y8vK6bJ+Lj3GlLJeTmpqqlJQUc7ugoEARERFXGhoAAACg1qlUYVFRUVHlJ77lllv02WefOe27//771alTJ82cOVMRERHy9vZWRkaGEhMTJf04ze2xY8cUGxsrSYqNjdWTTz6pEydOKDg4WJK0ceNGBQYGKjIy0uzz7rvvOp1n48aN5jF8fHwUHR2tjIwMs4CqqKhQRkaGJk2aJEmKjo6+YpbL8fX1la+vrzvDBAAAahiTQVQOk0Hgp37VAnlVISAgQF26dHHa17BhQzVt2tTcn5SUpJSUFDVp0kSBgYGaPHmyYmNjzVuw4uLiFBkZqZEjR2rhwoVyOBx69NFHNXHiRPML/fjx47VkyRLNmDFDY8eO1aZNm7Rq1SqlpaWZ501JSdHo0aN1/fXXq3fv3lq0aJGKiop0//33S5KCgoKumAUAAACoy35VYfHxxx/r/fff14kTJy65mvHXv/61SoJJ0rPPPiu73a7ExESVlJQoPj7e6fYrLy8vrV27VhMmTFBsbKwaNmyo0aNHO02F27ZtW6WlpSk5OVmLFy9Wy5Yt9fe//13x8fFmn6FDh+rkyZOaPXu2HA6HoqKilJ6e7vRA95WyAAAAAHWZy+tY/PnPf9ajjz6qjh07KiQkRDbb/7tcaLPZtGnTpioP6SlYxwKoPty6UDncugBcGZ8nlcPnSd3gyvdXl69YLF68WP/85z81ZsyYX5sPAAAAgIdxeYE8u92uG2+8sTqyAAAAAKilXC4skpOT9fzzz1dHFgAAAAC1lMu3Qj388MNKSEhQ+/btFRkZKW9vb6f2N998s8rCAQAAAKgdXC4spkyZovfff1/9+/dX06ZNnR7eBgAAAFA3uVxYLF++XG+88YYSEhKqIw8AAACAWsjlZyyaNGmi9u3bV0cWAAAAALWUy4XF3LlzNWfOHJ07d6468gAAAACohVy+Feq5557TV199pZCQELVp0+aSh7d3795dZeEAAAAA1A4uFxaDBg2qhhgAAAAAajOXC4s5c+ZURw4AAAAAtZjLz1gAAAAAwE+5fMXCbrf/4toV5eXlbgUCAAAAUPu4XFi89dZbTtvnz5/XJ598ouXLl+vxxx+vsmAAAAAAag+XC4u77777kn2DBw/Wddddp5UrVyopKalKggEAAACoParsGYs+ffooIyOjqg4HAAAAoBapksLihx9+0HPPPacWLVpUxeEAAAAA1DIu3wrVuHFjp4e3DcPQ2bNn1aBBA/373/+u0nAAAAAAageXC4tFixY5bdvtdjVv3lwxMTFq3LhxVeUCAAAAUIu4XFiMHj26OnIAAAAAqMVcLiwkKS8vTzt37tSJEydUUVHh1DZq1KgqCQYAAACg9nC5sHjnnXc0YsQIFRYWKjAw0Ol5C5vNRmEBAAAA1EEuzwr10EMPaezYsSosLFReXp7OnDljvk6fPl0dGQEAAABc5VwuLL777jtNmTJFDRo0qI48AAAAAGohlwuL+Ph47dq1qzqyAAAAAKilXH7GIiEhQdOnT9eBAwfUtWtXeXt7O7XfddddVRYOAAAAQO3gcmExbtw4SdK8efMuabPZbCovL3c/FQAAAIBaxeXC4qfTywIAAACAy89YAAAAAMBPUVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3/arCoqKiQp9//rm2bt2qDz74wOnlihdffFHdunVTYGCgAgMDFRsbq3Xr1pntxcXFmjhxopo2bSp/f38lJiYqNzfX6RjHjh1TQkKCGjRooODgYE2fPl1lZWVOfTZv3qyePXvK19dXHTp00LJlyy7J8vzzz6tNmzby8/NTTEyMdu7c6dRemSwAAABAXeVyYfHRRx+pQ4cO6ty5s2666Sb169fPfPXv39+lY7Vs2VILFixQVlaWdu3apZtvvll333239u/fL0lKTk7WO++8o9WrV2vLli3Kzs7WPffcY76/vLxcCQkJKi0t1fbt27V8+XItW7ZMs2fPNvscPXpUCQkJ6t+/v/bs2aNp06bpgQce0Pr1680+K1euVEpKiubMmaPdu3ere/fuio+P14kTJ8w+V8oCAAAA1GU2wzAMV94QFRWla6+9Vo8//rjCwsJks9mc2oOCgtwK1KRJEz399NMaPHiwmjdvrhUrVmjw4MGSpEOHDqlz587KzMxUnz59tG7dOt1xxx3Kzs5WSEiIJGnp0qWaOXOmTp48KR8fH82cOVNpaWnat2+feY5hw4YpLy9P6enpkqSYmBj16tVLS5YskfTjFZmIiAhNnjxZs2bNUn5+/hWzVEZBQYGCgoKUn5+vwMBAt8YJgDPb47Yrd4KMOS595AN1Ep8nlcPnSd3gyvdXl69YfPHFF/rzn/+szp07q1GjRgoKCnJ6/Vrl5eV6/fXXVVRUpNjYWGVlZen8+fMaMGCA2adTp05q1aqVMjMzJUmZmZnq2rWrWVRIUnx8vAoKCsyrHpmZmU7HuNDnwjFKS0uVlZXl1Mdut2vAgAFmn8pkuZySkhIVFBQ4vQAAAABP5HJhERMToy+//LLKAnz22Wfy9/eXr6+vxo8fr7feekuRkZFyOBzy8fFRo0aNnPqHhITI4XBIkhwOh1NRcaH9Qtsv9SkoKNAPP/yg77//XuXl5Zftc/ExrpTlcubPn+9UdEVERFRuUAAAAIBapp6rb5g8ebIeeughORwOde3aVd7e3k7t3bp1c+l4HTt21J49e5Sfn6//+7//0+jRo7VlyxZXY12VUlNTlZKSYm4XFBRQXAAAAMAjuVxYJCYmSpLGjh1r7rPZbDIMQzabTeXl5S4dz8fHRx06dJAkRUdH6+OPP9bixYs1dOhQlZaWKi8vz+lKQW5urkJDQyVJoaGhl8zedGGmpov7/HT2ptzcXAUGBqp+/fry8vKSl5fXZftcfIwrZbkcX19f+fr6ujAaAAAAQO3k8q1QR48eveR15MgR8z/dVVFRoZKSEkVHR8vb21sZGRlm2+HDh3Xs2DHFxsZKkmJjY/XZZ585zd60ceNGBQYGKjIy0uxz8TEu9LlwDB8fH0VHRzv1qaioUEZGhtmnMlkAAACAuszlKxatW7euspOnpqbqtttuU6tWrXT27FmtWLFCmzdv1vr16xUUFKSkpCSlpKSoSZMmCgwM1OTJkxUbG2vOwhQXF6fIyEiNHDlSCxculMPh0KOPPqqJEyeaVwrGjx+vJUuWaMaMGRo7dqw2bdqkVatWKS0tzcyRkpKi0aNH6/rrr1fv3r21aNEiFRUV6f7775ekSmUBAAAA6jKXCwtJ+uqrr7Ro0SIdPHhQkhQZGampU6eqffv2Lh3nxIkTGjVqlHJychQUFKRu3bpp/fr1GjhwoCTp2Wefld1uV2JiokpKShQfH68XXnjBfL+Xl5fWrl2rCRMmKDY2Vg0bNtTo0aM1b948s0/btm2Vlpam5ORkLV68WC1bttTf//53xcfHm32GDh2qkydPavbs2XI4HIqKilJ6errTA91XygIAQHVh+tPKYfpTwFour2Oxfv163XXXXYqKitKNN94oSdq2bZv27t2rd955xywKcCnWsQCqD1+8KocvXrUTf9+VU1V/34x35fB5Uje48v3V5SsWs2bNUnJyshYsWHDJ/pkzZ1JYAAAAAHWQyw9vHzx4UElJSZfsHzt2rA4cOFAloQAAAADULi4XFs2bN9eePXsu2b9nzx4FBwdXRSYAAAAAtYzLt0KNGzdODz74oI4cOaIbbrhB0o/PWDz11FNOi8EBAAAAqDtcLiwee+wxBQQE6JlnnlFqaqokKTw8XHPnztWUKVOqPCAAAACAq5/LhYXNZlNycrKSk5N19uxZSVJAQECVBwMAAABQe7j8jMXNN9+svLw8ST8WFBeKioKCAt18881VGg4AAABA7eByYbF582aVlpZesr+4uFgffvhhlYQCAAAAULtU+laoTz/91PznAwcOyOFwmNvl5eVKT09XixYtqjYdAAAAgFqh0oVFVFSUbDabbDbbZW95ql+/vv72t79VaTgAAAAAtUOlC4ujR4/KMAy1a9dOO3fuVPPmzc02Hx8fBQcHy8vLq1pCAgAAALi6VbqwaN26tSSpoqKi2sIAAAAAqJ1cnm5Wkr766istWrRIBw8elCRFRkZq6tSpat++fZWGAwAAAFA7uDwr1Pr16xUZGamdO3eqW7du6tatm3bs2KHrrrtOGzdurI6MAAAAAK5yLl+xmDVrlpKTk7VgwYJL9s+cOVMDBw6ssnAAAAAAageXr1gcPHhQSUlJl+wfO3asDhw4UCWhAAAAANQuLhcWzZs31549ey7Zv2fPHgUHB1dFJgAAAAC1jMu3Qo0bN04PPvigjhw5ohtuuEGStG3bNj311FNKSUmp8oAAAAAArn4uFxaPPfaYAgIC9Mwzzyg1NVWSFB4errlz52rKlClVHhAAAADA1c/lwsJmsyk5OVnJyck6e/asJCkgIKDKgwEAAACoPX7VOhYXUFAAAAAAkH7Fw9u5ubkaOXKkwsPDVa9ePXl5eTm9AAAAANQ9Ll+xGDNmjI4dO6bHHntMYWFhstls1ZELAAAAQC3icmGxdetWffjhh4qKiqqGOAAAAABqI5dvhYqIiJBhGNWRBQAAAEAt5XJhsWjRIs2aNUtff/11NcQBAAAAUBtV6laoxo0bOz1LUVRUpPbt26tBgwby9vZ26nv69OmqTQgAAADgqlepwmLRokXVHAMAAABAbVapwmL06NHVnQMAAABALVbpZywqKir01FNP6cYbb1SvXr00a9Ys/fDDD9WZDQAAAEAtUenC4sknn9Qjjzwif39/tWjRQosXL9bEiROrMxsAAACAWqLShcVrr72mF154QevXr9fbb7+td955R//5z39UUVFRnfkAAAAA1AKVLiyOHTum22+/3dweMGCAbDabsrOzqyUYAAAAgNqj0oVFWVmZ/Pz8nPZ5e3vr/Pnzv/rk8+fPV69evRQQEKDg4GANGjRIhw8fdupTXFysiRMnqmnTpvL391diYqJyc3Od+hw7dkwJCQlq0KCBgoODNX36dJWVlTn12bx5s3r27ClfX1916NBBy5YtuyTP888/rzZt2sjPz08xMTHauXOny1kAAACAuqhSs0JJkmEYGjNmjHx9fc19xcXFGj9+vBo2bGjue/PNNyt98i1btmjixInq1auXysrK9MgjjyguLk4HDhwwj5mcnKy0tDStXr1aQUFBmjRpku655x5t27ZNklReXq6EhASFhoZq+/btysnJ0ahRo+Tt7a0///nPkqSjR48qISFB48eP13/+8x9lZGTogQceUFhYmOLj4yVJK1euVEpKipYuXaqYmBgtWrRI8fHxOnz4sIKDgyuVBQAAAKirbIZhGJXpeP/991fqgK+++uqvDnPy5EkFBwdry5Ytuummm5Sfn6/mzZtrxYoVGjx4sCTp0KFD6ty5szIzM9WnTx+tW7dOd9xxh7KzsxUSEiJJWrp0qWbOnKmTJ0/Kx8dHM2fOVFpamvbt22eea9iwYcrLy1N6erokKSYmRr169dKSJUsk/TgLVkREhCZPnqxZs2ZVKsuVFBQUKCgoSPn5+QoMDPzV4wTgUrbHbVfuBBlzKvWRj6sMf9+VU1V/34x35fB5Uje48v210lcs3CkYKis/P1+S1KRJE0lSVlaWzp8/rwEDBph9OnXqpFatWplf5jMzM9W1a1ezqJCk+Ph4TZgwQfv371ePHj2UmZnpdIwLfaZNmyZJKi0tVVZWllJTU812u92uAQMGKDMzs9JZfqqkpEQlJSXmdkFBwa8dGgAAAOCqVulnLKpbRUWFpk2bphtvvFFdunSRJDkcDvn4+KhRo0ZOfUNCQuRwOMw+FxcVF9ovtP1Sn4KCAv3www/6/vvvVV5eftk+Fx/jSll+av78+QoKCjJfERERlRwNAAAAoHa5agqLiRMnat++fXr99detjlJlUlNTlZ+fb76OHz9udSQAAACgWlT6VqjqNGnSJK1du1YffPCBWrZsae4PDQ1VaWmp8vLynK4U5ObmKjQ01Ozz09mbLszUdHGfn87elJubq8DAQNWvX19eXl7y8vK6bJ+Lj3GlLD/l6+vr9LA7AAAA4KksvWJhGIYmTZqkt956S5s2bVLbtm2d2qOjo+Xt7a2MjAxz3+HDh3Xs2DHFxsZKkmJjY/XZZ5/pxIkTZp+NGzcqMDBQkZGRZp+Lj3Ghz4Vj+Pj4KDo62qlPRUWFMjIyzD6VyQIAAADUVZZesZg4caJWrFih//3vfwoICDCfVQgKClL9+vUVFBSkpKQkpaSkqEmTJgoMDNTkyZMVGxtrPiwdFxenyMhIjRw5UgsXLpTD4dCjjz6qiRMnmlcLxo8fryVLlmjGjBkaO3asNm3apFWrViktLc3MkpKSotGjR+v6669X7969tWjRIhUVFZmzYVUmCwAAAFBXWVpYvPjii5Kkfv36Oe1/9dVXNWbMGEnSs88+K7vdrsTERJWUlCg+Pl4vvPCC2dfLy0tr167VhAkTFBsbq4YNG2r06NGaN2+e2adt27ZKS0tTcnKyFi9erJYtW+rvf/+7uYaFJA0dOlQnT57U7Nmz5XA4FBUVpfT0dKcHuq+UBQAAAKirKr2OBdzHOhZA9WHe+cph3vnaib/vymEdi5rF50nd4Mr316tmVigAAAAAtReFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3WVpYfPDBB7rzzjsVHh4um82mt99+26ndMAzNnj1bYWFhql+/vgYMGKAvvvjCqc/p06c1YsQIBQYGqlGjRkpKSlJhYaFTn08//VS/+c1v5Ofnp4iICC1cuPCSLKtXr1anTp3k5+enrl276t1333U5CwAAAFBXWVpYFBUVqXv37nr++ecv275w4UI999xzWrp0qXbs2KGGDRsqPj5excXFZp8RI0Zo//792rhxo9auXasPPvhADz74oNleUFCguLg4tW7dWllZWXr66ac1d+5cvfzyy2af7du3695771VSUpI++eQTDRo0SIMGDdK+fftcygIAAADUVTbDMAyrQ0iSzWbTW2+9pUGDBkn68QpBeHi4HnroIT388MOSpPz8fIWEhGjZsmUaNmyYDh48qMjISH388ce6/vrrJUnp6em6/fbb9e233yo8PFwvvvii/vjHP8rhcMjHx0eSNGvWLL399ts6dOiQJGno0KEqKirS2rVrzTx9+vRRVFSUli5dWqkslVFQUKCgoCDl5+crMDCwSsYNwI9sj9usjlArGHOuio98uIi/78qpqr9vxrty+DypG1z5/nrVPmNx9OhRORwODRgwwNwXFBSkmJgYZWZmSpIyMzPVqFEjs6iQpAEDBshut2vHjh1mn5tuusksKiQpPj5ehw8f1pkzZ8w+F5/nQp8L56lMFgAAAKAuq2d1gJ/jcDgkSSEhIU77Q0JCzDaHw6Hg4GCn9nr16qlJkyZOfdq2bXvJMS60NW7cWA6H44rnuVKWyykpKVFJSYm5XVBQ8Av/jQEAAIDa66q9YuEJ5s+fr6CgIPMVERFhdSQAAACgWly1hUVoaKgkKTc312l/bm6u2RYaGqoTJ044tZeVlen06dNOfS53jIvP8XN9Lm6/UpbLSU1NVX5+vvk6fvz4Ff5bAwAAALXTVVtYtG3bVqGhocrIyDD3FRQUaMeOHYqNjZUkxcbGKi8vT1lZWWafTZs2qaKiQjExMWafDz74QOfPnzf7bNy4UR07dlTjxo3NPhef50KfC+epTJbL8fX1VWBgoNMLAAAA8ESWFhaFhYXas2eP9uzZI+nHh6T37NmjY8eOyWazadq0aXriiSe0Zs0affbZZxo1apTCw8PNmaM6d+6sW2+9VePGjdPOnTu1bds2TZo0ScOGDVN4eLgkafjw4fLx8VFSUpL279+vlStXavHixUpJSTFzTJ06Venp6XrmmWd06NAhzZ07V7t27dKkSZMkqVJZAAAAgLrM0oe3d+3apf79+5vbF77sjx49WsuWLdOMGTNUVFSkBx98UHl5eerbt6/S09Pl5+dnvuc///mPJk2apFtuuUV2u12JiYl67rnnzPagoCBt2LBBEydOVHR0tJo1a6bZs2c7rXVxww03aMWKFXr00Uf1yCOP6JprrtHbb7+tLl26mH0qkwUAAACoq66adSzqAtaxAKoP885XDvPO1078fVcO61jULD5P6gaPWMcCAAAAQO1BYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxGYQEAAADAbRQWAAAAANxWz+oAqDm2x21WR6gVjDmG1REAAABqHa5YAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt1FYAAAAAHAbhQUAAAAAt9WzOgAAAADqLtvjNqsj1ArGHMPqCFfEFQsXPf/882rTpo38/PwUExOjnTt3Wh0JAAAAsBxXLFywcuVKpaSkaOnSpYqJidGiRYsUHx+vw4cPKzg42Op4AFBj+IWxcmrDL4wAUFW4YuGCv/71rxo3bpzuv/9+RUZGaunSpWrQoIH++c9/Wh0NAAAAsBRXLCqptLRUWVlZSk1NNffZ7XYNGDBAmZmZl31PSUmJSkpKzO38/HxJUkFBQfWG/TnF1py2tqmq/32C5gdVyXE8XX5qftUciL/vSqmyzx/Gu1IY75rFeNcsxrtmWfX98cJ5DePKV2BtRmV6QdnZ2WrRooW2b9+u2NhYc/+MGTO0ZcsW7dix45L3zJ07V48//nhNxgQAAACq3PHjx9WyZctf7MMVi2qUmpqqlJQUc7uiokKnT59W06ZNZbNxf3JBQYEiIiJ0/PhxBQYGWh3H4zHeNYvxrlmMd81ivGsW412zGG9nhmHo7NmzCg8Pv2JfCotKatasmby8vJSbm+u0Pzc3V6GhoZd9j6+vr3x9fZ32NWrUqLoi1lqBgYH8H7cGMd41i/GuWYx3zWK8axbjXbMY7/8nKKhyt3fz8HYl+fj4KDo6WhkZGea+iooKZWRkON0aBQAAANRFXLFwQUpKikaPHq3rr79evXv31qJFi1RUVKT777/f6mgAAACApSgsXDB06FCdPHlSs2fPlsPhUFRUlNLT0xUSEmJ1tFrJ19dXc+bMueR2MVQPxrtmMd41i/GuWYx3zWK8axbj/esxKxQAAAAAt/GMBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBuFBQAAAAC3UVgAAAAAcBvrWKDGnDp1Sp9++qm6d++uJk2a6Pvvv9c//vEPlZSUaMiQIercubPVEQEAAPArsY4FasTOnTsVFxengoICNWrUSBs3btSQIUNUr149VVRUKDs7W1u3blXPnj2tjupxNm3apK1btyonJ0d2u13t2rXTXXfdpWuuucbqaB6pvLxcXl5e5vaOHTtUUlKi2NhYeXt7W5jMs7zxxhu67bbb1KBBA6uj1FlHjx7Vl19+qbCwMHXp0sXqOB6nqKhIWVlZTp/dPXv2lM1mszqaRzpy5Mgl/64cOHCgAgMDrY5WuxhADRgwYIDxwAMPGAUFBcbTTz9ttGzZ0njggQfM9vvvv98YNGiQhQk9T25urtG7d2/Dbrcb9erVM+x2uxEdHW2EhoYaXl5exvTp062O6FGys7ONG2+80fDy8jJuuukm4/Tp00ZCQoJhs9kMm81mXHvttUZ2drbVMT2GzWYzAgMDjXHjxhkfffSR1XE83oQJE4yzZ88ahmEY586dMxITEw273W7YbDbDbrcb/fv3N9vhnvLycmP69OlGgwYNDLvdbo6zzWYzWrdubaxZs8bqiB6lsLDQGDx4sDnGdrvd/Pekv7+/sWTJEqsj1io8Y4EakZWVpZSUFAUEBGjq1KnKzs7WuHHjzPZJkybp448/tjCh55kyZYrCw8N15swZFRYW6g9/+IOuu+465eTkaMOGDfrnP/+pxYsXWx3TY8ycOVOGYeitt95SWFiY7rjjDhUUFOj48eP6+uuv1bx5cz355JNWx/QoDz/8sHbt2qXY2Fh16dJFixYt0qlTp6yO5ZFeeuklnTt3TpL0pz/9STt27NB7772nwsJCffDBBzp27Bh/31XkkUce0dq1a7Vy5UqtX79effv21YIFC3TgwAGNGjVKQ4YM0YYNG6yO6TFSUlKUk5OjTz/9VJ9//rnuuecejRo1SgUFBVq8eLFmzJihFStWWB2z9rC6skHd0LBhQ+Po0aPmtr+/v/HVV1+Z2998843h5+dnQTLPFRgYaOzbt8/cLiwsNLy9vY38/HzDMAzjX//6l9GxY0er4nmcsLAwIzMz0zAMwzh16pRhs9mM9957z2zPyMgw2rVrZ1U8j2Oz2Yzc3FzDMAxj165dxoQJE4xGjRoZvr6+xpAhQ4wNGzZYnNCzXDzeXbp0MVasWOHU/r///c+49tprrYjmccLCwowPPvjA3P72228Nf39/o7i42DAMw5g3b54RGxtrVTyP06xZM2PXrl3m9unTpw0/Pz+jqKjIMAzDWLJkiREVFWVVvFqHKxaoERERETpy5Ii5/frrryssLMzczsnJUbNmzayI5rF8fX2d7sW12+0qLy9XWVmZJOmGG27Q119/bVE6z3PmzBm1aNFCktSkSRM1aNBArVu3Nts7dOignJwcq+J5tOjoaL3wwgvKycnRK6+8opMnT+rWW29V27ZtrY7mUS58njgcDnXr1s2prXv37jp+/LgVsTxOYWGh+VkiSWFhYSouLtaZM2ckSYmJidq7d69V8TxOWVmZ03MU/v7+KisrU1FRkSQpLi5Ohw4dsiperUNhgRoxbNgwnThxwtxOSEhQ/fr1ze01a9aod+/eVkTzWH379tXs2bNVVFSk8+fP65FHHlG7du3UpEkTSdLJkyfVuHFji1N6juDgYKfCYdKkSeZYSz8WHg0bNrQimke63AOsfn5+GjlypN5//30dPnxYw4cPtyCZ53rssceUkpIiu92u7Oxsp7ZTp07x911Funbtqv/+97/m9qpVq+Tv76/Q0FBJUkVFhXx9fa2K53F69erldFvw4sWL1bx5czVv3lzSj4Wev7+/VfFqHaabRY2YM2fOL7b/8Y9/dJpJB+77y1/+ori4ODVq1Eg2m00NGzbU6tWrzfaDBw9qzJgx1gX0MFFRUcrMzDQL5AULFji1b9269ZJfefHrGVeY0LBDhw7c81+FbrrpJh0+fFiSFBkZqW+++cap/d1339V1111nRTSPM2/ePCUkJGjNmjXy8/PT9u3b9fTTT5vt6enp6tGjh4UJPcuCBQs0cOBAvfHGG/Lx8ZHD4dDy5cvN9u3bt+v222+3MGHtwnSzgAc7d+6ctm3bppKSEvXp04fbzSy0c+dONWjQgGk5q8g333yjVq1aMfXmVeLIkSPy8fFRy5YtrY7iEfbu3atVq1appKRE8fHxGjhwoNWRPFpOTo7Wrl2rkpIS3XzzzYqMjLQ6Uq1FYYGrwv/+9z/l5+dr1KhRVkcBAADAr0BhgatCp06d9MUXX6i8vNzqKHXGrl27dO7cOd10001WR6kTcnJydP78ebVq1crqKHUCf981i/GGp+Kz2zUUFkAd1blzZ33++ecUczWE8a5ZjHfNYrxrDmNdsxhv1/DwNlBHZWRk6Pz581bHqDNee+01c4ExVD/+vmsW411z5s+fr/z8fKtj1Bl8druGKxawRF5enlavXq1jx46pdevWGjJkiIKCgqyOBQAAgF+JwgI14p577tHw4cM1ePBg7d+/X/369ZPNZlO7du309ddfy2azadOmTercubPVUT2Ow+HQjh075HA4JEmhoaGKiYkx50RH9Tp//ry8vb2tjuGxioqKlJWVpZycHNntdrVr1049e/ZktqhqVlZWpvfff9/8cah///5MGV7Nli1bpt/97nf8CFeNjh075vRZ0rRpU6sj1T6WrPeNOqdx48bGwYMHDcMwjNtuu80YPny4UVJSYhiGYZSWlhpJSUlGXFyclRE9TmFhoTFixAjDy8vLqFevnhEcHGwEBwcb9erVM7y8vIz77rvPKCoqsjqmx1i5cqX5N20YhvG3v/3NaNWqlWG3242mTZsajz/+uIXpPE95ebkxffp0o0GDBobdbjfsdrths9kMm81mtG7d2lizZo3VET3KpEmTjHfeeccwDMM4fvy40alTJ8PLy8sICQkxvLy8jK5duxrffvutxSk9m7e3t3HgwAGrY3ik559/3vy8vvh14403Grt27bI6Xq3CytuoEcXFxeavtnv27NHDDz8sHx8fSZK3t7dmzJihHTt2WBnR40ydOlU7d+5UWlqaiouLlZubq9zcXBUXF+vdd9/Vzp07NXXqVKtjeox7771XeXl5kqRXX31V06dP15gxY/TOO+8oOTlZCxcu1N///ndrQ3qQRx55RGvXrtXKlSu1fv169e3bVwsWLNCBAwc0atQoDRkyRBs2bLA6psdYvXq12rRpI0l66KGH1LJlSzkcDjkcDp04cUKtW7fWtGnTLM3oKZo0aXLZV1lZmWJjY81tVI2//OUvevLJJzV9+nS99NJL6tixo+bOnau0tDS1a9dON910k3bt2mV1zFqDW6FQI/r06aOkpCSNGzdOPXv21OzZszVo0CCzfePGjRo1apRycnKsC+lhGjdurLS0NN1www2Xbd+2bZvuuOMOnTlzpoaTeSa73S6Hw6Hg4GDFxMRo8ODBmj59utn+4osv6pVXXtHu3bstTOk5wsPDtXLlSv3mN7+RJH333Xfq1KmTvv/+e/n6+upPf/qT1q1bp+3bt1uc1DPUr19fBw4cUNu2bRUREaE33njDXGVekvbt26f+/fvr5MmTFqb0DAEBAfrtb3+rIUOGmPsMw9ADDzygefPmqUWLFpKk0aNHWxXRo7Rt21YvvPCCbrvtNknS559/rhtuuEEOh0P16tXT1KlTdfDgQX6oqCRmhUKNeOyxxzRq1Ch5e3trypQpSk5O1qlTp9S5c2cdPnxYc+bM0ciRI62O6VEqKirMq0KX4+Pjo4qKihpM5Pku3Nd/5MgRxcXFObXFxcVp5syZVsTySIWFheYXLEkKCwtTcXGxzpw5o9DQUCUmJmrBggUWJvQs1157rXbu3Km2bdsqICBABQUFTu1nz57l86SKfPLJJxo+fLg2bdqk559/Xv7+/pKkcePGadCgQawKXcVOnDjh9HznNddco/z8fJ08eVJhYWEaO3as+vbta2HC2oVboVAjEhIS9PLLL+uxxx5TUlKSvvnmG40bN059+/bVH/7wByUmJmr+/PlWx/Qod9xxhx588EF98sknl7R98sknmjBhgu68804Lknmu9PR0rVmzRn5+fpdMT1hcXMwDxVWoa9eu+u9//2tur1q1Sv7+/uakBBUVFfL19bUqnsdJTk7Www8/rM2bNys1NVVTpkxRRkaGsrOz9f777+v3v/+97rnnHqtjeoQOHTpo+/btCg0NVVRUlLZt22Z1JI927bXXauPGjeb2+++/Lx8fH/OzxM/Pj89uF3DFAjUmMTFRgwYNUlZWlo4ePaqKigqFhYUpOjpaAQEBVsfzOEuWLNHw4cMVHR2txo0bKzg4WNKPv87k5eUpPj5eS5YssTilZ7n41oRNmzYpNjbW3P7oo4/Uvn17K2J5pHnz5ikhIcEs5LZv366nn37abE9PT1ePHj0sTOhZxowZo9OnTyshIUGGYai8vNzpqtxdd92lZ5991sKEnqVevXp66qmnFB8fr+HDh2vEiBF8ua0mqampuu+++/Tee+/Jz89Pb775pqZMmWKO9+bNm9WlSxeLU9YePGMBeLiDBw/qo48+cppuNjY2Vp06dbI4Wd2ydu1aeXt7Kz4+3uooHmPv3r1atWqVSkpKFB8fr4EDB1odyePl5eVp48aNOnLkiPnj0I033qhrrrnG6mge69SpUxo3bpzef/99ffTRR+rYsaPVkTzOunXr9O9//9v8LBk3bpzZdurUKUli6tlKorBAjdq0aZO2bt3qNE/0XXfdxb+UAAAAajkKC9SIEydO6M4779SuXbtkt9tVUVGhHj166LvvvtPJkyeVkpKihQsXWh3TI1HMXR0uLOR20003WR3Foxw5cuSSv++BAwcqMDDQ6mgeLS8vT6tXrzYXyBsyZAgLt1WTixcjbNOmjfr168dihLh6WbWABuqWoUOHGoMGDTLy8/ON4uJiY9KkScaoUaMMwzCMjIwMo2nTpsaiRYssTulZcnNzjd69ext2u92oV6+eYbfbjejoaCM0NNTw8vIypk+fbnXEOmXPnj2G3W63OobHKCwsNAYPHmwuime3282/bX9/f2PJkiVWR/Qov/vd74zVq1cbhmEY+/btM5o1a2Y0b97ciImJMUJCQozQ0FAWb6siLEZYs0pLS43p06cb7du3N3r16mX84x//cGp3OBx8druAWaFQI9atW6cnnnhCgYGB8vX11YIFC/Tf//5XBQUFuvnmm7Vo0SK9+OKLVsf0KFOmTFF4eLjOnDmjwsJC/eEPf9B1112nnJwcbdiwQf/85z+1ePFiq2MCv0pKSopycnL06aef6vPPP9c999yjUaNGqaCgQIsXL9aMGTO0YsUKq2N6jIsfYJ0+fbri4uL07bff6qOPPtLx48eVkJDAAnlVhMUIa9aTTz6p1157TePHj1dcXJxSUlL0+9//3qmPwc09lWd1ZYO6oXnz5sb+/fvN7XPnzhl2u904deqUYRiG8dVXXxm+vr5WxfNIgYGBxr59+8ztwsJCw9vb28jPzzcMwzD+9a9/GR07drQqnsdp3LjxL74CAwP51asKNWvWzNi1a5e5ffr0acPPz88oKioyDMMwlixZYkRFRVkVz+PUr1/f+PLLLw3DMIywsDBj9+7dTu2HDx82goKCLEjmefz8/IwjR44YhmEYLVu2NHbs2OHU/tlnnxnNmjWzIppH6tChg3mFyDAM44svvjA6dOhgjBkzxqioqOCKhYuYbhY1om/fvpo9e7aWL18uHx8fPfLII2rXrp2aNGkiSTp58qQaN25scUrP4uvr6zQ9od1uV3l5ucrKyiRJN9xwg77++muL0nmekpISTZgwQV27dr1s+zfffKPHH3+8hlN5rrKyMqfnKPz9/VVWVqaioiI1aNBAcXFxevjhhy1M6Fm6deumTZs2qX379goNDdU333zjNJ3vN998o/r161uY0HOwGGHN+u6775ymk+3QoYM2b96sm2++WSNHjuT5TxdRWKBG/OUvf1FcXJwaNWokm82mhg0bavXq1Wb7wYMHNWbMGOsCeiCKuZoVFRWliIgIp7UsLrZ3714KiyrUq1cvLV682FyLZfHixWrevLmaN28u6ceVuS+sWAz3PfbYYxo1apS8vb01ZcoUJScn69SpU+rcubMOHz6sOXPmaOTIkVbH9AgXFiMMCQkxFyP829/+Zo711KlTWYywCoWGhuqrr74ybz+TpBYtWuj9999X//79+W7iIgoL1Ih27drp008/1datW1VaWqo+ffqoWbNmMgxDNpuN/+NWA4q5mpWQkKC8vLyfbW/SpIlGjRpVc4E83IIFCzRw4EC98cYb8vHxkcPh0PLly8327du36/bbb7cwoWdJSEjQyy+/rGnTpik7O1uGYZhz/fv6+mr8+PGaP3++xSk9A4sR1qybb75ZK1as0C233OK0Pzw8XJs2bVK/fv2sCVZLMd0sLOXj46O9e/eqc+fOVkfxSOfOndO2bdtUUlJiFnOAp8jJydHatWtVUlKim2++WZGRkVZH8njl5eXavXu30wJ50dHRCggIsDqax8nLy9OGDRt09OhRFiOsRt98840OHTr0s4uXZmdna+PGjT97NRrOKCxQI1JSUi67f/HixbrvvvvMFS3/+te/1mQsAAAAVBFuhUKNWLRokbp3765GjRo57TcMQwcPHlTDhg2dHjRG9cvNzdVLL72k2bNnWx3Fo3z77bdq1KjRJff3nz9/XpmZmSyQV4VOnTqlTz/9VN27d1eTJk30/fff6x//+IdKSko0ZMgQroRWI8MwtHnzZn355ZcKCwtTfHy8vL29rY7lUXbu3KnMzEw5HA5JPz4LEBsbq969e1uczLOUlJTIbrebf79fffWV/vnPf5qLPyYlJalt27YWp6xFLJuPCnXK/PnzjbZt2xoZGRlO++vVq+c0DS1qDgu2Va3s7GyjV69eht1uN7y8vIyRI0caZ8+eNduZsrBq7dixwwgKCjJsNpvRuHFjY9euXUbbtm2Na665xmjfvr1Rv359Iysry+qYHuO2224z8vLyDMMwjFOnThkxMTGGzWYzmjdvbtjtdqNTp07GiRMnLE7pGXJzc42+ffsaNpvNaN26tdG7d2+jd+/eRuvWrQ2bzWb07dvXyM3NtTqmx/jtb39rLv64detWw9fX1+jWrZsxdOhQo0ePHkaDBg2M7du3W5yy9uBWKNSYjz/+WPfdd5/uvPNOzZ8/X97e3vL29tbevXu5N7oafPrpp7/YfujQId17770qLy+voUSebfTo0Tp8+LCWLFmivLw8zZo1SzabTRs2bFDjxo2Vm5ursLAwpomsIgMHDlSbNm3017/+VS+99JIWL16sW2+9Va+88ookaezYsTpz5ozeeusti5N6BrvdLofDoeDgYP3hD3/Qli1btHbtWrVt21bffvutBg0apF69erHQaRUYPHiwsrOz9eqrr6pjx45ObYcPH9bYsWMVHh7uNBkHfr2goCDt2rVL11xzjfr166eePXs63Zb92GOP6f3339fWrVstTFmLWF3ZoG45e/asMWrUKKNbt27GZ599Znh7e3PFoprYbDbDbrcbNpvtkteF/fyCXnXCw8OdFrIqLi427rzzTiMqKso4deoUVyyqWOPGjY0DBw4YhmEYpaWlht1udxr/rKwso0WLFlbF8zg2m838lbxjx47G//73P6f29957z2jbtq0V0TyOv7//JQsQXmzXrl2Gv79/DSbybA0bNjQOHjxoGIZhhISEGHv27HFq//LLLxlvF9itLmxQt/j7+2v58uVKTU3VgAED+LW8GjVp0kSvvPKKjh49esnryJEjWrt2rdURPUp+fr7TuiC+vr5688031aZNG/Xv318nTpywMJ3nKS0tNRdk8/b2VoMGDZxmPWvWrJlOnTplVTyPdOE5uDNnzqh9+/ZObR06dFB2drYVsTyOr6/vJYviXezs2bPy9fWtwUSeLSYmRu+8844kqX379tq7d69T+549e8z1n3BlPLwNSwwbNkx9+/ZVVlaWWrdubXUcjxQdHa3s7OyfHd+8vDwZ3AlZZS6s1XLxVJD16tXT6tWrNWTIEN1xxx0WpvM8EREROnLkiLmo1euvv66wsDCzPScnh+mVq9iYMWPk6+ur8+fP6+jRo7ruuuvMNofDccnkHPh1hg4dqtGjR+vZZ5/VLbfcYq4wX1BQoIyMDKWkpOjee++1OKXneOKJJ3TbbbepqKhI9957rx566CF98cUX5oKEzz33nFJTU62OWWtQWMAyLVu2VMuWLa2O4bHGjx+voqKin21v1aqVXn311RpM5Nluu+02vfzyy0pMTHTaf6G4SExM1PHjxy1K53mGDRvmdBUoISHBqX3NmjXMnlOFLp7D/+6779a5c+ec2t944w1FRUXVcCrP9Ne//lUVFRUaNmyYysrK5OPjI+nHq3T16tVTUlKS/vKXv1ic0nPExsZq3bp1SklJ0Y4dOyRJTz75pKQfF8mbO3eupk6damXEWoWHtwGgCpSVlencuXPmr4uXa//uu++4QldDzp07Jy8vL24ZqSFFRUXy8vKSn5+f1VE8RkFBgbKyspymm42Ojv7Zzxi47+TJk06LP164IorK4xkLoI46fvy4xo4da3UMj1GvXr1f/Bd+Tk6OHn/88RpMVLedOnVKEyZMsDpGnXH69Gn94Q9/sDqGxzh48KDeeOMNhYWF6d5771WPHj20atUqTZs2TZs2bbI6nsc5ePCgXn31VZ0+fVoxMTFq3LixnnrqKY0dO5bxdhFXLIA6au/everZsycP0NcQxrtmMd41i/GuOunp6br77rvl7++vc+fO6a233tKoUaPUvXt3VVRUaMuWLdqwYYNuvvlmq6N6BMa7avGMBeCh1qxZ84vtR44cqaEkdQPjXbMY75rFeNecefPmafr06XriiSf0+uuva/jw4ZowYYJ5339qaqoWLFjAF90qwnhXLa5YAB7KbrfLZrP94sxPNpuNXxirCONdsxjvmsV415ygoCBlZWWpQ4cOqqiokK+vr3bu3KkePXpIkvbt26cBAwaYz17APYx31eIZC8BDhYWF6c0331RFRcVlX7t377Y6okdhvGsW412zGO+adWHNELvdLj8/PwUFBZltAQEBys/PtyqaR2K8qw6FBeChoqOjlZWV9bPtV/r1Ea5hvGsW412zGO+a06ZNG33xxRfmdmZmplq1amVuHzt2zGnNFriH8a5aPGMBeKjp06f/4joWHTp00Pvvv1+DiTwb412zGO+axXjXnAkTJjjdUtalSxen9nXr1nG/fxVivKsWz1gAAAAAcBu3QgEAAABwG4UFAAAAALdRWAAAAABwG4UFAAAAALdRWAAAAABwG4UFAABXmd/97ndq3LixBg8ebHUUAKg0CgsAAK4yU6dO1WuvvWZ1DABwCYUFAAD/v1OnTik4OFhff/21pTn69eungICAS/YPGzZMzzzzjAWJAODKKCwAwAONGTNGNpvtkteXX35pdbSr2pNPPqm7775bbdq0sTrKZT366KN68sknlZ+fb3UUALhEPasDAACqx6233qpXX33VaV/z5s0v6VdaWiofH5+ainXVOnfunP7xj39o/fr11X6uqKgolZWVXbJ/w4YNCg8P/9n3denSRe3bt9e///1vTZw4sTojAoDLuGIBAB7K19dXoaGhTi8vLy/169dPkyZN0rRp09SsWTPFx8dLkioqKjR//ny1bdtW9evXV/fu3fV///d/TscsKirSqFGj5O/vr7CwMD3zzDPq16+fpk2bZvZp06aNFi1a5PS+qKgozZ07t9Ln6devn6ZMmaIZM2aoSZMmCg0NNd9/QUVFhRYuXKgOHTrI19dXrVq10pNPPqnXXntNTZs2VUlJiVP/QYMGaeTIkT87Xu+++658fX3Vp0+fK57j4pyTJ0/WtGnT1LhxY4WEhOiVV15RUVGR7r//fgUEBKhDhw5at26d07n27Nmjffv2XfL6paLigjvvvFOvv/76FfsBQE2jsACAOmj58uXy8fHRtm3btHTpUknS/Pnz9dprr2np0qXav3+/kpOTdd9992nLli3m+6ZPn64tW7bof//7nzZs2KDNmzdr9+7dLp27Mue5kLFhw4basWOHFi5cqHnz5mnjxo1me2pqqhYsWKDHHntMBw4c0IoVKxQSEqIhQ4aovLxca9asMfueOHFCaWlpGjt27M/m+vDDDxUdHe207+fO8dOczZo1086dOzV58mRNmDBBQ4YM0Q033KDdu3crLi5OI0eO1Llz51wap5/Tu3dv7dy585LCCQAsZwAAPM7o0aMNLy8vo2HDhuZr8ODBhmEYxm9/+1ujR48eTv2Li4uNBg0aGNu3b3fan5SUZNx7772GYRjG2bNnDR8fH2PVqlVm+6lTp4z69esbU6dONfe1bt3aePbZZ52O0717d2POnDmVOs+FjH379nXq06tXL2PmzJmGYRhGQUGB4evra7zyyiuX/e8/YcIE47bbbjO3n3nmGaNdu3ZGRUXFZfsbhmHcfffdxtixY83tK53jcjnLysqMhg0bGiNHjjT35eTkGJKMzMzMnz3OT91yyy1Gs2bNjPr16xstWrRwGq+9e/cakoyvv/660scDgJrAMxYA4KH69++vF1980dxu2LCh+c8//WX+yy+/1Llz5zRw4ECn/aWlperRo4ck6auvvlJpaaliYmLM9iZNmqhjx46VzlSZ81zQrVs3p+2wsDCdOHFCknTw4EGVlJTolltuuex5xo0bp169eum7775TixYttGzZMvOB9p/zww8/yM/Pz9y+0jkul9PLy0tNmzZV165dzX0XrnBcyF4Z77333s+21a9fX5Kq7AoIAFQVCgsA8FANGzZUhw4dfrbtYoWFhZKktLQ0tWjRwqnN19fXpfPa7XYZhuG07/z58y6fx9vb22nbZrOpoqJC0v/7cv1zevTooe7du+u1115TXFyc9u/fr7S0tF98T7NmzXTmzBlz+0rn+KWcF++7UMxcyO6u06dPS7r8g/gAYCWesQAAKDIyUr6+vjp27Jg6dOjg9IqIiJAktW/fXt7e3tqxY4f5vjNnzujzzz93Olbz5s2Vk5NjbhcUFOjo0aOVPk9lXHPNNapfv74yMjJ+ts8DDzygZcuW6dVXX9WAAQOuePwePXrowIEDLp3DCvv27VPLli3VrFkzq6MAgBOuWAAAFBAQoIcffljJycmqqKhQ3759lZ+fr23btikwMFCjR4+Wv7+/kpKSNH36dDVt2lTBwcH64x//KLvd+Teqm2++WcuWLdOdd96pRo0aafbs2fLy8qr0eSrDz89PM2fO1IwZM+Tj46Mbb7xRJ0+e1P79+5WUlCRJGj58uB5++GG98sorlVrFOj4+XqmpqTpz5owaN25cqXNY4cMPP1RcXJxl5weAn0NhAQCQJP3pT39S8+bNNX/+fB05ckSNGjVSz5499cgjj5h9nn76aRUWFurOO+9UQECAHnrooUsWa0tNTdXRo0d1xx13KCgoSH/605/MKxaVPU9lPPbYY6pXr55mz56t7OxshYWFafz48WZ7UFCQEhMTlZaWpkGDBl3xeF27dlXPnj21atUq/f73v6/UOWpacXGx3n77baWnp1uWAQB+js346Y2wAAC4oF+/foqKirpk7YqrwS233KLrrrtOzz33XKX6p6Wlafr06dq3b98lV2KuBi+++KLeeustbdiwweooAHAJrlgAADzOmTNntHnzZm3evFkvvPBCpd+XkJCgL774Qt99951Lz3zUFG9vb/3tb3+zOgYAXBaFBQDA4/To0UNnzpzRU0895dJ0uJKcVhG/2jzwwANWRwCAn8WtUAAAAADcdvXdQAoAAACg1qGwAAAAAOA2CgsAAAAAbqOwAAAAAOA2CgsAAAAAbqOwAAAAAOA2CgsAAAAAbqOwAAAAAOA2CgsAAAAAbqOwAAAAAOA2CgsAAAAAbqOwAAAAAOA2CgsAAAAAbvv/AAlIiqRvDj/nAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample4 = []\n",
    "for ele in sample3.items():\n",
    "    # print(ele[0].state)\n",
    "    sample4.append(ele[0].state * ele[1])\n",
    "counts4 = np.sum(sample4, axis=0)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.ylabel(\"Photon number\")\n",
    "plt.xlabel(r\"Frequency (cm$^{-1}$)\")\n",
    "plt.xticks(range(len(omegap)), np.round(omegap, 1), rotation=90)\n",
    "plt.bar(range(len(omegap)), counts4, color=\"green\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "piquasso",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
